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Even  simple,  deterministic  rules  can  generate  interesting  behavior  in  dynamical 
systems.  This  dissertation  examines  some  real  world  systems  for  which  fairly  simple, 
locally  defined  rules  yield  useful  or  interesting  properties  in  the  system  as  a  whole.  In 
particular,  we  study  routing  in  peer-to-peer  networks  and  the  motion  of  crystal  steps. 

Peers  can  vary  by  three  orders  of  magnitude  in  their  capacities  to  process  net¬ 
work  traffic.  This  heterogeneity  inspires  our  use  of  “proportionate  load  balancing,” 
where  each  peer  provides  resources  in  proportion  to  its  individual  capacity.  We  pro¬ 
vide  an  implementation  that  employs  small,  local  adjustments  to  bring  the  entire 
network  into  a  global  balance.  Analytically  and  through  simulations,  we  demon¬ 
strate  the  effectiveness  of  proportionate  load  balancing  on  two  routing  methods  for 
de  Bruijn  graphs,  introducing  a  new  “reversed”  routing  method  which  performs  better 
than  standard  forward  routing  in  some  cases. 

The  prevalence  of  peer-to-peer  applications  prompts  companies  to  locate  the 


hosts  participating  in  these  networks.  We  explore  the  use  of  supervised  machine 
learning  to  identify  peer-to-peer  hosts,  without  using  application-specific  informa¬ 
tion.  We  introduce  a  model  for  “triples,”  which  exploits  information  about  nearly 
contemporaneous  flows  to  give  a  statistical  picture  of  a  host’s  activities.  We  find  that 
triples,  together  with  measurements  of  inbound  vs.  outbound  traffic,  can  capture 
most  of  the  behavior  of  peer-to-peer  hosts. 

An  understanding  of  crystal  surface  evolution  is  important  for  the  development 
of  modern  nanoscale  electronic  devices.  The  most  commonly  studied  surface  features 
are  steps,  which  form  at  low  temperatures  when  the  crystal  is  cut  close  to  a  plane 
of  symmetry.  Step  bunching,  when  steps  arrange  into  widely  separated  clusters  of 
tightly  packed  steps,  is  one  important  step  phenomenon.  We  analyze  a  discrete  model 
for  crystal  steps,  in  which  the  motion  of  each  step  depends  on  the  two  steps  on  either 
side  of  it.  We  find  an  time-dependence  term  for  the  motion  that  does  not  appear  in 
continuum  models,  and  we  determine  an  explicit  dependence  on  step  number. 
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Chapter  1 


Introduction 

1.1  Overview  of  the  dissertation 

Even  simple,  deterministic  rules  can  generate  interesting  behavior  in  dynamical  sys¬ 
tems.  For  example,  the  logistic  map  exhibits  chaos  for  certain  parameter  values.  This 
dissertation  examines  some  real-world  systems  for  which  fairly  simple,  locally-defined 
rules  yield  useful  or  interesting  properties  in  the  system  as  a  whole. 

The  remaining  sections  of  this  chapter  provide  an  introduction  to  modern  peer- 
to-peer  networks  and  to  crystal  surfaces.  We  trace  the  history  of  peer-to-peer  net¬ 
works  from  completely  centralized  systems  like  Napster,  to  completely  decentralized 
systems  like  Gnutella,  to  hybrid  systems  like  Fast  Track,  to  alternative  systems  like 
Bit  Torrent.  We  then  give  a  detailed  introduction  to  distributed  hash  tables,  which  are 
decentralized  but  structured  peer-to-peer  systems.  Finally,  we  discuss  the  motivation 
behind  the  study  of  crystal  steps  and  the  step  bunching  phenomenon. 

In  Chapter  2,  we  address  the  issue  of  heterogeneity  in  peer-to-peer  networks. 
By  heterogeneity,  we  refer  to  the  fact  that  peers  can  vary  by  three  orders  of  mag¬ 
nitude  in  their  capacities  to  process  traffic.  We  define  and  analyze  networks  based 
on  two  routing  methods  for  de  Bruijn  graphs,  introducing  a  new  “reversed”  routing 
method  which  performs  better  than  standard  forward  routing  in  some  cases.  We 
quantify  how  heterogeneity  affects  the  average  query  path  length  under  each  routing 
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method,  proving  that  the  length  decreases  in  most  cases.  We  also  provide  an  imple¬ 
mentation  of  what  we  term  “proportionate  load  balancing,”  which  is  applicable  to 
a  range  of  peer-to-peer  networks.  Proportionate  load  balancing  employs  small  local 
adjustments  to  bring  the  entire  network  into  a  global  balance,  where  each  participant 
provides  resources  in  proportion  to  its  capacity.  We  demonstrate  the  effectiveness  of 
proportionate  load  balancing  on  de  Bruijn  networks. 

Chapter  3  explores  the  use  of  supervised  machine  learning  to  identify  peer- 
to-peer  hosts,  without  using  port  numbers  or  other  application-specific  information. 
Most  studies  in  this  area  have  attempted  to  classify  TCP  flows,  rather  than  individual 
users.  We  introduce  a  model  for  “triples,”  which  capture  the  composite  behavior 
available  when  we  study  hosts  rather  than  flows.  Triples  exploit  information  about 
nearly-contemporaneous  flows  to  give  a  statistical  picture  of  a  host  as  a  whole.  We 
find  that  triples,  together  with  measurements  of  inbound  vs.  outbound  traffic,  can 
capture  most  of  the  behavior  of  peer-to-peer  hosts,  while  maintaining  acceptable 
“false-positive”  rates. 

In  Chapter  4,  we  analyze  a  discrete  model  for  crystal  steps.  In  this  model,  the 
motion  of  each  step  depends  on  the  two  steps  on  either  side  of  it.  We  evaluate  the 
effect  of  step-step  interactions  on  the  stability  of  the  crystal  surface,  with  a  focus  on 
step  bunching. 
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1.2  Peer-to-peer  networks 


Classical  internet  activity  fits  the  client-server  model,  where  powerful  central  comput¬ 
ers  (the  servers)  provide  resources  to  lower-capacity  end  users  (the  clients).  Some  ex¬ 
amples  of  servers  include  cnn.com,  which  disseminates  news;  microsoft.com,  which 
distributes  Windows  updates;  yahoo  .  com,  which  provides  forums  for  games  and  chat¬ 
ting;  or  the  machines  that  handle  all  umd .  edu  email. 

Peer-to-peer  (P2P)  applications,  on  the  other  hand,  represent  a  distributed  type 
of  internet  activity.  A  peer-to-peer  network  is  a  collection  of  machines  that  associate 
with  one  another  to  share  files  or  other  resources.  A  P2P  network  may  have  servers 
that  facilitate  these  communications,  but  unlike  in  a  client-server  setup,  the  servers 
do  not  provide  all  of  the  content. 

1.2.1  Statistics 

Up  to  60%  of  worldwide  Internet  traffic  is  peer-to-peer  [42],  Most  peer-to-peer  ap¬ 
plications  are  used  for  file  sharing  —  sending  files  directly  between  end-users.  By 
volume,  P2P  traffic  is  [42] 

•  61%  video  content, 

•  11%  audio  content, 

•  11%  compressed  files, 

•  17%  other. 
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(a)  A  centralized  P2P  system. 


(b)  A  decentralized  P2P  system. 


Figure  1.1:  Peer-to-peer  applications  may  be  divided  broadly  into  two  classes. 
Centralized  systems  have  a  server  that  organizes  file  sharing  for  a  large  number 
of  subsidiary  users;  this  server  is  a  bottleneck  and  single  point  of  failure  for 
the  system.  Decentralized  systems  lack  a  central  authority,  so  they  avoid  the 
weaknesses  of  centralized  systems;  however,  searches  can  only  follow  a  limited 
number  of  paths,  so  such  systems  can  be  highly  inefficient. 

$2.3  billion  in  copyrighted  movies  were  downloaded  worldwide  in  2005  [36].  Only  13% 
of  that  loss  came  from  users  in  the  United  States. 


1.2.2  History 

Broadly,  there  are  two  classes  of  peer-to-peer  applications:  centralized  and  decentral¬ 
ized.  See  Figure  1.1. 

Centralized  systems  (Napster,  DirectConnect,  BitTorrent)  have  a  computer  that 
acts  as  a  hub  to  organize  hie  transmissions  for  a  large  number  of  client  computers. 
Napster  had  a  single  hub  for  all  its  users;  DirectConnect  has  many  disjoint  sets  of 
clients  and  a  hub;  BitTorrent  has  a  hub  for  every  individual  hie.  Despite  their  central- 
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ized  organization,  these  systems  are  peer-to-peer  because  shared  hies  are  transmitted 
directly  between  clients.  The  hubs  act  only  as  facilitators  to  help  clients  find  other 
users  who  possess  the  information  they  want,  not  as  content  servers.  The  hubs  also 
provide  an  obvious  target  for  disabling  the  network;  Napster  was  shut  down  by  court 
order  in  July  2001  [53]. 

Decentralized  systems  (Gnutella,  Fast  Track  (Kazaa),  distributed  hash  tables 
(Kademilia))  have  each  peer  connected  to  some  small  number  of  other  peers,  without 
a  central  authority,  so  searches  for  hies  must  be  distributed  through  their  networks. 
Gnutella  became  popular  with  the  demise  of  Napster.  It  employs  hooding,  where  a 
peer  asks  its  connected  peers  if  they  have  a  certain  hie,  whereupon  they  ask  their 
connected  peers,  etc.  Fast  Track  automatically  divides  its  peers  into  clients  and  su¬ 
pernodes,  with  supernodes  performing  searches  in  a  Gnutella-like  fashion  for  their 
connected  clients.  Distributed  hash  tables  are  structured  to  provide  more  efficient 
searches  in  exhange  for  the  overhead  of  maintaining  the  structure.  In  all  cases,  once 
the  requested  hie  is  located,  it  is  transferred  between  the  relevant  clients. 

Some  systems  (eDonkey2000)  seem  to  lie  between  the  two  extremes.  eDon- 
key2000  originally  acted  like  DirectConnect  with  disjoint  hub  networks,  but  it  now 
allows  the  hubs  to  pass  searches  among  themselves. 
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1.3  Distributed  hash  tables 


Distributed  hash  tables  (DHTs)  form  completely  decentralized  P2P  systems  which, 
due  to  an  imposed  structure,  have  reliable  and  efficient  search  methods.  See  [32]  for 
details  on  several  popular  protocols.  Let  V  be  the  set  of  peers.  In  general,  DHTs 
have  three  structural  components: 

1.  A  keyspace  /C,  with  a  publicly-known  function  /  :  {data  items}  — >  /C  that 
assigns  keys  to  data  items. 

2.  A  rule  for  assigning  a  portion,  or  zone,  of  the  keyspace  to  each  peer  i  £  V,  such 
that  1J  ieVzone{i)  =  /C. 

3.  A  set  of  application-level  connections,  assigned  between  peers  based  on  the  keys 
in  their  zones. 

In  the  following  subsections,  we  formalize  this  definition,  then  briefiy  discuss 
the  operation  and  maintenance  of  DHTs. 

1.3.1  Continuous-discrete  definition 

For  definiteness,  we  describe  Naor  and  Wieder’s  continuous-discrete  approach  [37]  to 
defining  DHTs.  In  this  framework, 

1.  /C  is  a  continuous  space.  There  is  a  graph  Gcont  that  has  K.  as  its  vertex  set.  (/ 
is  unrestricted;  it  is  often  taken  to  be  a  uniform  hash  function.) 

2.  /C  is  decomposed  into  possibly-overlapping,  connected  zones,  one  per  peer. 
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3.  If  GCOnt  has  an  edge  from  zone(i )  to  zone(j),  then  there  may  be  an  application- 
level  link  from  peer  i  to  peer  j .  This  yields  a  discrete  graph  GdiSC  on  the  vertex 
set  V  of  peers. 

In  Chapter  2,  we  will  define  proportionate  load  balancing  for  DHTs  that  have 
a  one-dimensional  keyspace  /C  =  [0, 1)  mod  1.  With  this  type  of  DHT,  peer  zones 
are  simply  intervals.  Peer  i  has  label  id(i),  and  should  keep  track  of  (at  least)  its 
immediate  predecessor  pred(i)  and  successor  succ(i),  which  are  the  peers  such  that 

id(j )  G  (. id(pred(i)),id(succ(i )))  j  =  i.  (1.1) 

We  use  the  following  definitions  in  this  paper. 

•  Zone:  zone(i)  =  [id(i) ,  id(succ(i)))  modi 

•  Zone  size:  Z{  =  id(succ(i ))  —  id(i )  mod  1 

•  Relative  zone  size:  rt  =  nzi,  the  ratio  of  zt  to  the  average  zone  size  ^ 

1.3.2  Data  sharing 

DHTs  operate  successfully  because  every  peer  knows  to  treat  the  location  idx  = 
f(x)  G  /C  as  the  authority  on  information  about  the  data  item  x.  Suppose  idx  e 
zone(i).  Fundamentally,  DHTs  need  to  provide  just  one  operation:  lookup  :  /C  — >  V  : 
idx  i — >  i,  which  passes  messages  through  the  network  to  locate  the  peer  in  charge  of 
item  x. 

Data  items  (hies)  enter  the  network  by  having  a  peer  elect  to  share  them.  If 
peer  o  (the  “owner”)  elects  to  share  item  x,  it  executes  lookup(f(x))  to  find  the  peer 
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i  (the  “intermediary”)  in  charge  of  idx.  Peer  o  then  gives  its  contact  information  to 
peer  i. 

When  another  peer  s  (the  “seeker”)  wants  to  find  item  x,  it  also  executes 
lookup(f  (x))  to  find  the  peer  i.  Peer  i  reports  peer  o’s  contact  information  to  peer  s, 
so  that  s  can  obtain  x  from  o.  If  no  peer  has  elected  to  share  the  item  a;,  peer  i  can 
definitively  alert  peer  s  to  this  fact. 

In  a  network  of  n  peers,  many  (though  not  all)  DHTs  provide  search  times  of 
order  O  (log2(n))  with  O  (log2(n))  connections  per  peer.  In  Chapter  2,  we  will  explain 
the  lookup  operation  in  detail  for  de  Bruijn  networks.  These  networks  provide  search 
times  of  order  O  (log2(n))  with  only  O  ( k )  connections  per  peer. 

1.3.3  Maintenance  of  the  network 

New  peers  join  the  network  by  a  process  known  as  bootstrapping.  A  new  peer  must 
locate  a  peer  already  in  the  network,  often  via  a  list  of  IP  addresses  for  computers 
known  to  have  been  connected  in  the  past.  The  new  peer  locates  a  target  through 
this  intermediary,  and  consults  with  the  target  to  divide  that  zone  between  them. 

A  peer  can  sign  off  of  the  network  by  transferring  its  zone  to  another  peer  with 
a  consecutive  region  of  the  keyspace.  If,  instead,  the  peer  simply  drops  out  of  the 
network  without  notifying  any  of  its  neighbors,  a  period  of  time  is  necessary  for  the 
network  to  completely  recover.  Ensuring  satisfactory  performance  during  this  period 
is  a  subject  of  much  study  [44].  Typically,  peers  will  re-advertise  their  available  hies 
periodically  in  order  to  assist  in  the  recovery. 


1.4  Crystal  surfaces 


An  understanding  of  crystal  surface  evolution  is  important  for  modern  nanoscale 
applications,  such  as  nanowires  and  quantum  dots  [69].  A  variety  of  surface  features 
[81]  can  occur.  For  sufficiently  low  temperatures,  crystal  steps  are  the  dominant 
feature  on  a  crystal  surface  cut  close  to  a  plane  of  symmetry. 

Specifically,  a  crystal  lattice  has  certain  planes  of  symmetry  that  depend  on  the 
lattice  structure.  A  cubic  lattice  has  three  perpendicular  planes  of  symmetry.  The 
crystal  can  be  cut  at  a  “vicinal  angle”  (f>  «  1  away  from  a  plane  of  symmetry.  At 
high  temperatures,  the  surface  will  be  statistically  rough,  but  at  temperatures  T  less 
than  the  “roughening  temperature”  TR  (which  depends  on  the  material),  observable 
steps  form.  See  Figure  1.2.  The  terraces  have  an  average  length  L  oc  cot (</>). 

Step  bunching  is  a  widely  studied  surface  phenomenon.  In  the  one-dimensional 
model,  one  equilibrium  state  is  a  train  of  uniformly-spaced  steps.  In  contrast,  under 
the  step  bunching  instability,  steps  arrange  themselves  into  widely-separated  clusters 
of  tightly-packed  steps.  This  bunching  has  been  observed  in  a  variety  of  experimental 
systems,  most  commonly  as  a  result  of  electromigration  due  to  an  applied  direct 
current  [70],  but  also  as  a  result  of  material  deposition  onto  the  surface  [65].  Such 
experimental  observations  have  motivated  interest  in  a  theoretical  understanding  of 
the  phenomenon. 
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(a)  Crystal  steps  with  curvature. 


(b)  Straight  crystal  steps. 


Figure  1.2:  At  sufficiently  low  temperatures,  steps  are  observable  on  a  crystal 
surface  cut  close  to  a  plane  of  symmetry.  The  steps  are  approximately  one 
lattice  constant  high,  and  the  terraces  are  typically  tens  or  hundreds  of  lattice 
constants  wide.  Real  crystal  steps  have  kinks,  islands,  vacancies,  and  other 
imperfections,  but  they  can  be  effectively  modeled  by  smooth  curved  steps  as 
in  Figure  1.2(a).  Under  certain  conditions,  the  steps  can  be  considered  to  have 
zero  curvature,  as  in  Figure  1.2(b).  This  allows  one-dimensional  modeling  of 


step  positions. 


Chapter  2 


Exploiting  heterogeneity 

2.1  Introduction 

Internet-enabled  computers  are  heterogeneous,  in  that  their  capacities  to  process 
traffic  vary  widely  [46,  29].  In  a  peer-to-peer  situation,  it  can  be  considered  equitable 
for  the  more  capable  peers  to  handle  a  larger  fraction  of  the  traffic  [44],  We  will 
equalize  the  utilizations  of  the  peers,  where  utilization  is  defined  in  Section  2.3.1  as 
the  fraction  of  a  peer’s  capacity  that  is  used  by  the  P2P  application. 

We  refer  to  the  operation  of  equalizing  utilizations  as  proportionate  load  balanc¬ 
ing  (PLB).  This  name  distinguishes  it  from  typical  notions  of  load  balancing,  where 
the  absolute  amount  of  traffic  to  each  peer  is  equalized.  We  will  demonstrate  analyt¬ 
ically  and  through  simulations  that  leveraging  heterogeneity  can  improve  the  query 
path  length  and  congestion  in  a  particular  class  of  distributed  hash  tables,  the  de 
Bruijn  DHTs. 

2.1.1  Related  work 

DHTs  based  on  the  de  Bruijn  architecture  have  been  studied  previously  [16,  37,  25, 
31,  17,  2],  We  will  present  the  basics  of  both  forward  and  reversed  de  Bruijn  net¬ 
works.  In  particular,  we  offer  an  algorithm  for  routing  in  reversed  de  Bruijn  networks 
(“reversed  routing”)  that  is  novel  for  being  provably  correct,  i.e.,  the  destination  is 
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always  reached  in  absence  of  node  failures.  Broose  [17]  provided  a  reversed  routing 
algorithm  where  the  probability  of  failing  to  reach  the  desired  destination  was  nonzero 
(but  small  under  the  assumption  of  homogeneous  zone  sizes).  Distance  Halving  [37] 
defined  its  underlying  architecture  in  a  reversed  manner,  but  only  provided  routing 
algorithms  that  relied  (in  whole  or  in  part)  on  forward  routing.  D2B  [16],  Koorde 
[25],  and  ODRI  [31]  worked  exclusively  with  forward  networks.  All  of  these  employed 
an  assumption  of  homogeneity  in  the  partitioning  of  the  keyspace. 

Demand  on  the  network’s  peers  can  differ  for  several  reasons:  (i)  nonuniform 
numbers  of  data  items  may  be  assigned  to  the  peers  [21];  (ii)  data  items  of  varying 
popularities  may  receive  highly  nonuniform  numbers  of  requests  [28];  and  (iii)  the  net¬ 
work  structure  may  cause  demand  to  vary  inherently,  as  we  will  see  in  Section  2. 3. 2. 2. 
Some  prior  approaches  to  coping  with  peer  heterogeneity  have  divided  the  peers  into 
a  hierarchy  of  two  or  more  levels  based  on  their  capacities  [19].  Other  approaches 
have  employed  a  varying  number  of  “virtual  servers”  [9,  22]  for  each  peer.  These 
approaches  essentially  use  (i)  to  address  heterogeneous  peer  capacities,  but  lack  the 
real-time  adaptability  to  deal  with  (ii)  and  (iii). 

Our  PLB  implementation  allows  continuous  variability  in  load  distribution,  and 
keeps  the  simplicity  of  a  single-level  DHT  structure.  Unlike  a  recent  implementation 
[20],  it  allows  continuous  variability  in  peer  capacities,  and  seeks  to  balance  the  ratio 
of  these  capacities  to  the  actual  traffic  load  at  each  peer,  not  simply  to  the  amount 
of  keyspace  assigned  to  each  peer.  It  also  operates  throughout  the  peers’  lifetimes, 
not  just  during  the  process  of  joining  or  leaving  the  network. 

Initially,  proportionate  load  balancing  seems  similar  to  load  balancing  in  het- 
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erogeneous  processor  networks  [14].  While  our  implementation  of  PLB  is  similar  to 
first-order  diffusion  schemes  for  processor  networks,  or  the  NbrAdjust  operation  for 
parallel  databases  [18],  the  problem  setting  is  different.  In  P2P  networks,  traffic  load 
is  intimately  tied  to  network  structure.  We  modify  the  structure  to  balance  queries 
at  all  stages  along  the  search  path,  not  just  in  the  final  stage  of  actually  serving  the 
response  to  a  request. 

2.1.2  DHTs  and  proportionate  load  balancing 

As  noted  above,  most  protocols  call  for  all  peers  to  hold  zones  of  approximately  equal 
size  [54],  so  that  absolute  traffic  levels  are  nearly  equal.  However,  by  assigning  more 
traffic  to  peers  with  higher  capacity,  PLB  can  balance  their  relative  traffic  levels. 
We  will  demonstrate  how  this  is  possible  in  certain  DHTs,  where  the  number  of 
connections  a  peer  has,  and  thus  the  amount  of  traffic  it  processes,  depends  on  the 
size  of  its  zone. 

All  DHTs  that  use  the  continuous-discrete  framework  of  Section  1.3.1  satisfy 
the  following  assumption,  which  is  critical  for  our  implementation  of  proportionate 
load  balancing. 

(Al)  The  expected  number  of  routing  paths  that  pass  through  a  given  peer  in¬ 
creases  with  the  relative  size  of  its  zone. 

Well-known  examples  where  Assumption  Al  holds  include  CAN  [32]  and  ODRI  [31]. 
The  proofs  in  Section  2.2.2  will  rely  in  part  on  the  following  assumption. 

(A2)  The  peers  occur  in  a  random  order  in  the  keyspace,  that  is,  P(j  =  succ(i))  = 
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— l-r  for  all  j  ^  i. 

n—  1  J  ' 

2.1.3  Overview 

The  remainder  of  this  chapter  is  structured  as  follows.  In  Section  2.2  we  define  two 
classes  of  de  Brnijn  DHTs  and  analytically  examine  them  under  heterogeneous  struc¬ 
ture.  In  Section  2.3  we  present  a  simple  scheme  for  dynamically  adapting  DHTs  (not 
necessarily  de  Brnijn)  to  exploit  heterogeneity.  In  Section  2.4  we  present  simulation 
results. 


14 


2.2  Heterogeneous  de  Bruijn  networks 


Researchers  have  studied  de  Bruijn  DHTs  because  they  offer  the  potential  for  loga¬ 
rithmic  diameter  logfc (Speers)  with  constant  degree  k.  We  use  them  because  they 
have  node  degree  intimately  tied  to  zone  size  [54],  making  them  good  candidates  for 
proportionate  load  balancing. 

Some  arguments  have  been  made  against  the  practicality  of  de  Bruijn  DHTs  [10]. 
First,  heterogeneous  partitioning  of  the  keyspace  prevents  de  Bruijn  peers  from  having 
constant  degree  [10].  We  embrace  this  disparity  by  providing  the  larger  zones  to  peers 
with  more  capacity.  Second,  some  peers  in  a  DHT  must  have  degree  f2(log(#peers)) 
(unlike  the  peers  in  a  classical  de  Bruijn  network)  if  the  network  is  to  remain  connected 
when  half  the  peers  fail  [25].  By  allowing  varying  peer  degrees  (Section  2.2.2. 1)  and 
using  an  average  degree  k  >  2  ([10]  footnote  1),  we  can  achieve  a  balance  between 
fault-tolerance  and  maintenance  costs.  Third,  the  path  between  two  peers  may  be 
fl(^peers)  in  the  worst  case.  If  the  ordering  of  the  peers  in  the  keyspace  is  random, 
the  probability  of  this  is  low;  furthermore,  heterogeneity  improves  the  expected  query 
path  length  (Section  2. 2. 2. 2). 

In  this  section,  we  will  present  the  basics  of  both  forward  and  reversed  de  Bruijn 
networks.  The  authors  of  Koorde  and  D2B  provided  more  complete  P2P  protocols 
[25,  16]  than  we  offer  here.  The  authors  of  ODRI  offered  comparisons  of  de  Bruijn  with 
other  DHT  structures,  including  discussions  of  diameter,  routing,  and  fault  tolerance 

[31]- 
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2.2.1  Forward  and  reversed  de  Bruijn  networks 

In  this  section,  we  first  define  continuous  versions  of  forward  and  reversed  de  Bruijn 
graphs.  We  then  show  how  to  discretize  these  systems  to  obtain  de  Bruijn  distributed 
hash  tables. 

2. 2. 1.1  de  Bruijn  graphs  (continuous  model) 

The  underlying  space  is  the  unit  interval  [0,1)  mod  1.  When  considering  graphs  of 
degree  k  (integer  >  2),  we  will  represent  the  values  in  base  k,  i.e., 

OO 

a  =  .010203  •  •  •  =  —  e  [0, 1)  (2.1) 

i=  1 

where  each  a*  G  Z*.  =  (0, 1,  * . . ,  k  —  1}.  The  value  of  k  should  be  a  reasonable  number 
of  connections  for  each  peer  to  maintain,  e.g.,  k  =  8  or  k  =  16. 

2. 2. 1.1.1  Forward  de  Bruijn  graphs  Algebraically,  each  point  a  G  [0,1)  has 
one  outgoing  edge,  to  the  point  ka  mod  1.  Symbolically,  this  is  represented  as  an 
edge  .0x0203  •  •  •  — >  .020304  •  •  • .  We  can  thus  refer  to  these  as  left-shiftmg  graphs. 

2. 2. 1.1. 2  Reversed  de  Bruijn  graphs  Algebraically,  each  point  a  G  [0,1)  has 
k  outgoing  edges,  to  the  points  (a  +  x)/k  mod  1,  for  x  G  Z*,.  Symbolically,  these 
are  represented  as  edges  ,aia2a3  ■  ■  ■  — >  .xaia2  ■  ■  ■  ■  We  can  thus  refer  to  these  as 
right- shifting  graphs. 

2. 2. 1.1. 3  Comments  The  outgoing  neighbors  of  a  vertex  in  a  reversed  de  Bruijn 
graph  are  its  incoming  neighbors  in  the  corresponding  forward  de  Bruijn  graph,  and 
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.010 


.001 


.01 - ^.00 


.10 - *.11 


(a)  k  =  D  =  2 


(b)  k  =  2,  D  =  3 


Figure  2.1:  Two  reversed  de  Bruijn  graphs  on  finite  spaces  {.0102  ■  ■  ■  od  :  a*  G 
Zfc}.  Each  graph  has  kD  vertices;  each  vertex  has  in-  and  out-degree  k  =  2. 

The  diameter  of  each  graph  is  observably  D. 

vice  versa. 

In  some  applications,  the  underlying  space  may  be  taken  as  a  collection  of 
elements  {.ci^  -  ■  ■  ar>  :  a*  G  Zfc}  (D  fixed)  that  is  finite  but  still  much  larger  (e.g., 
2128)  than  the  set  of  peers  might  ever  be.  In  this  case,  each  vertex  in  the  forward 
graph  has  k  outgoing  edges:  .0102  •  •  •  au-i — >  .0203  •  •  •  aoy,  y  G  Z&.  The  diameter 
(the  maximum  length  of  any  shortest  path  between  vertices)  of  the  resultant  forward 
and  reversed  graphs  is  D.  Fig.  2.1  shows  two  reversed  de  Bruijn  graphs  on  small 
underlying  spaces. 


2.2. 1.2  de  Bruijn  networks  (discrete  model) 

As  in  Section  1.3,  to  create  a  distributed  hash  table,  we  assign  each  peer  i  an  iden¬ 
tifier  id(i )  in  accordance  with  Eq.  (1.1).  Peer  i  then  controls  the  segment  zone(i )  = 
[. id(i ),  id(succ(i))),  and  has  links  to  all  peers  in  control  of  zones  where  some  underlying 


17 


.01010 


.01010 


.01111 


.00001 


.10010 


.10110 


.11010 


.10110 


.11010 


(a)  Outgoing  edges 


(b)  Incoming  edges 


Figure  2.2:  Here  k  =  2.  The  six  large  dots  are  the  ids  chosen  by  six  peers, 
and  the  arcs  depict  the  zone  each  peer  is  in  charge  of.  Fig.  2.2(a)  shows 
the  outgoing  edges  for  the  peer  with  key  .010100,  while  Fig.  2.2(b)  shows 
the  incoming  edges  for  that  same  peer.  The  dashed  lines  are  some  of  the 
relevant  edges  in  the  underlying  reversed  de  Bruijn  graph.  The  peer  will  have 
network-level  connections  to  (or  from)  all  zones  matching  at  least  one  of  these 
de  Bruijn  edges  —  the  solid  lines  show  these  resulting  connections  in  the  P2P 
network. 

edge  from  zone(i)  terminates.  See  Fig.  2.2. 

We  now  describe  algorithms  for  routing  in  de  Bruijn  DHTs,  and  verify  their 
correctness.  We  will  suppose  that  peer  i  initiates  a  request  for  an  item  with  key 
b  =  .&1&2&3  •  •  ■ ,  which  lies  in  some  other  peer’s  zone(j). 

2.2. 1.2.1  Forward  routing  Let  d  be  such  that  k~d  <  Zi.  Then  for  any  sequence 
.61&2&3  •  •  • ,  there  exist  aq,  a2, . . . ,  G  7Lk  such  that  .Gqaq  •  •  •  ddb^bz  •  •  •  G  zone(i). 
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At  most  two  values  of  (ai,  a2, . . . ,  ad)  are  required,  and  peer  i  can  calculate  them. 
Following  the  sequence 

.aia2'  •  •  ad- 1  ad  b\  b2  ■  ■  — > 

.a2U3' •  •  ad  b\  62  b3  ■  ■  — : ► 

(2.2) 

,adb\  -  ■  ■  bd-2bd-i  bd  bd+i - > 

• bib2 ■  ■  ■  bd- 1  bd  bd+1bd+2 •  ■  ■  =  6 

of  underlying  edges  yields  a  valid  path  of  peer  links  through  the  network,  starting 
at  peer  i  and  terminating  at  the  appropriate  peer  j  in  d  steps.  The  next  hop  is 
computable  locally  by  each  peer  involved. 

2. 2. 1.2. 2  Reversed  routing  Let  d  be  such  that  ( b  —  k~d,  b  +  k~d )  C  zone(j).1 
(Peer  %  cannot  determine  d  under  our  assumption  of  heterogeneity;  we  will  deal  with 
this  difficulty  soon.)  Then  any  point  .bib2  ■  ■  ■  bdCiC2c3  ■  ■  ■  lies  in  ( b  —  k~d ,  b  +  k~d)  and 
thus  is  covered  by  zone(j ) .  Following  any  sequence 

■  Cl C2 •  '  '  Cd  Cd+lCd+2'  - - * 

■bdCl'  ‘  ‘  Crf— 1  Cd  Cd+ 1'  '  '  * 

(2.3) 

,b2b3-  •  •  ci  c2  c3 - > 

. M2 •  •  •  bd  Ci  c2  ■  ■  ■  e  zone(j) 

1What  happens  if  b  =  id(j )?  In  the  continuous  case,  this  happens  with  probability  0,  but  can 
still  be  corrected  by  requiring  peer  j  —  1  to  forward  requests  for  b  to  peer  j.  In  the  case  of  a  finite 
underlying  space  {.aia2  •  •  •  am  :  a*  €  Z*,},  we  may  take  d  =  D. 
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of  underlying  edges  yields  a  valid  path  of  peer  links,  starting  at  an  arbitrary  point  c 
and  terminating  at  the  appropriate  peer  j  in  d  steps. 

Since  the  optimal  value  of  d,  dm in,  cannot  be  determined  a  priori  by  the  source 
peer  i,  we  must  modify  the  procedure  in  Eq.  (2.3).  Peer  i  should  perform  exponential 
polling  until  d,m-m  is  exceeded,  trying  values  d  =  1,2,4,8,16,....  The  queries  need 
not  be  returned  to  i  between  failed  attempts.  That  is,  using  to  denote  the  d-step 
procedure  in  Eq.  (2.3),  and  starting  at  id(i)  =  a  =  .a3a2a3  •  •  • ,  the  sequence 


a 

.b\a 

.bib2b\d 

.bib2b3bj)ib2bid 


,b\  a 
.bib2  bia 
■bifobsbi  bib2bid 
.&!•••  b8  bib2b3b4bib2bid 


(2.4) 


reaches  zone(j)  in  finite  time.  Letting  2h~1  <  dmin  <  2h,  the  total  number  of  steps 
is£  =  l  +  2  +  4  +  --  -  +  2h  =  2h+1  —  1  G  [2dmin  —  1, 4dmin  —  1).  The  next  hop  is 
computable  locally  by  each  peer  involved,  as  long  as  the  current  “stage”  (1,  2,  4,  8, 
. . . )  is  included  with  the  message  header. 

Alternatively,  peer  i  can  send  out  multiple  messages  simultaneously,  a  .bid, 
d  ^ 2  .bib2d,  . . . ,  a  ^>29  -b\  ■  ■  ■  b2gd.  If  no  response  is  received,  g  should  be  increased. 
The  value  of  g  can  be  based  on  prior  experience.  The  number  of  messages  is  the  same 
as  Eq.  (2.4),  but  the  response  time  is  shorter  because  some  are  sent  in  parallel. 

Of  course,  in  either  case,  it  is  probably  desirable  to  have  a  larger  first  attempt 

d  b2fd,  2^  >  1. 
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2.2.2  Heterogeneity  in  de  Bruijn  networks 


We  now  offer  analytical  analyses  of  de  Bruijn  DHTs  under  heterogeneity.  To  our 
knowledge,  this  is  the  first  such  analysis  for  either  forward  or  reversed  networks. 
Section  2.2.2. 1  shows  that  forward  networks  are  much  more  liable  to  develop  peers 
with  only  one  outgoing  neighbor,  which  are  susceptible  to  disconnection  from  the 
graph.  Section  2. 2. 2. 2  demonstrates  that  high  variance  among  zone  sizes  always 
reduces  the  query  path  length  of  reversed  networks  under  reasonable  assumptions,  but 
that  there  are  reasonable  assumptions  for  forward  networks  that  yield  longer  query 
path  lengths.  Finally,  Section  2. 2. 2. 3  discusses  the  use  of  caching  as  a  complementary 
tool  to  PLB. 

2. 2. 2.1  Number  of  neighbors 

We  begin  with  a  simple  property:  node  degree.  If  we  were  assuming  homogeneity 
and  roughly  equal  zone  sizes,  all  peers  would  have  in-degree  and  out-degree  very  near 
k  [31].  However,  as  the  zone  sizes  vary  in  our  heterogeneous  networks,  the  degrees 
become  non-uniform. 

2. 2. 2. 1.1  Forward  and  reversed  neighbors 

Theorem  2.1.  Suppose  peer  i  has  relative  zone  size  rt  in  either  a  forward  or  reversed 
de  Bruijn  DHT  with  n  peers. 

•  Averaging  over  forward  de  Bruijn  DHTs,  peer  i  expects  to  have  k{wd  =  1  +  kri 
outgoing  links  to  other  peers. 
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•  Averaging  over  reversed  de  Bruijn  DHTs,  peer  i  expects  to  have  kfev  =  k  +  r8 

outgoing  links  to  other  peers. 

Note  that  k{wd  arid  k™  are  independent  of  n. 

Proof.  We  begin  with  a  general  fact:  If  a  DHT  has  n  arbitrarily-ordered  peers,  and 
[a,  b )  is  an  interval  in  the  keyspace  [0, 1),  then  the  expected  number  of  distinct  peers 
in  charge  of  some  or  all  of  [a,  b)  is  n{b  —  a)  +  1.  Since  the  n  peers  are  arranged  in 
no  particular  order  by  Assumption  A2,  the  expected  number  of  peers  that  lie  inside 
(a,  b)  is  n(b  —  a).  Each  of  these  peers  is  in  charge  of  a  portion  of  [a,  b).  Additionally, 
there  is  another  peer  whose  key  either  is  a  or  immediately  precedes  a;  this  peer  is  in 
charge  of  the  first  portion  of  [a,  b). 

Now  recall  that  zone(i )  =  [id(i),  id(succ(i ))). 

In  the  forward  case,  zone[i )  maps  to  [kid(i),  kid(succ(i))),  an  interval  of  length 
kzi .  Then  peer  i  expects  to  connect  to  n  ■  kzi  +  1  =  1  +  krt  peers. 

In  the  reversed  case,  zone(i )  maps  to  k  evenly-spaced  intervals  ]d(su^cM) )  _|_ 
|,  x  G  Zfc.  Each  of  these  intervals  has  length  Then  peer  i  expects  to  connect  to 
k(n  ■  fzi  +  1)  =  k  +  Ti  peers.  □ 

In  the  reversed  case,  even  if  is  very  small,  krtev  ~  k.  Therefore,  even  a  very 
weak  peer  will  have  alternative  paths  to  route  a  query  through  if  one  of  its  neighbors 
fails.  In  the  forward  case,  if  r*  is  very  small,  then  /qfwd  ps  1  and  the  peer  will  be 
susceptible  to  disconnection  from  the  graph  ([54],  Section  2). 

The  expected  number  of  links  for  a  peer  with  an  average-sized  zone  (r*  =  1)  is 
k  +  1,  rather  than  k,  which  is  the  number  of  links  in  a  uniform  de  Bruijn  graph.  This 
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is  because  the  uniform  graph  is  unstable  with  respect  to  all  of  its  zone  boundaries 
lining  up. 

We  note  that  the  number  of  incoming  links  in  a  forward  DHT  is  the  number  of 
outgoing  links  in  the  corresponding  reversed  DHT,  and  vice  versa.  However,  incoming 
links  do  not  directly  help  with  routing  around  a  failed  neighbor. 

2. 2. 2. 1.2  Cost  of  PLB  In  Section  2.3.3,  we  will  discuss  the  number  of  peers  who 
must  be  notified  when  peer  i  participates  in  our  implementation  of  PLB.  In  the  worst 
case,  all  the  in-  and  out-neighbors  of  either  peer  pred(i)  or  i  must  be  notified  of  a 
change  in  id(i). 

By  Theorem  2.1,  we  expect  the  number  of  neighbors  that  must  be  contacted  (in 
either  the  forward  or  reversed  case)  to  be  <  kiwd  +  kTev  =  (k+  l)(max{rj,  rpre(j(i)}  + 1). 

2. 2. 2. 2  Query  path  lengths 

If  we  were  assuming  homogeneity,  so  all  n  peers  had  roughly  equal  zone  sizes,  then 
we  would  have  d  ~  logfc  n  for  both  the  forward  and  reversed  cases.  Thus  the  number 
of  routing  steps  would  be  C(logfcn),  as  in  [37,  25,  16,  17,  31].  We  note  that  linear 

polling  in  the  reversed  case  (trying  values  d  =  1,  2,  3,4,  5, . . .)  would  have  resulted  in 

0{{\ogkn)2)  steps. 

In  this  paper  we  assume  heterogeneity  of  peers  and  zone  sizes.  The  average 
number  of  routing  steps  will  depend  on  these  quantities: 

•  r  =  (ri, _ ,  rn),  the  relative  zone  sizes. 

•  p  —  (pi, _ ,  pn ),  the  relative  rates  at  which  the  peers  request  items.  Reasonable 
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values:  p  =  1  (uniform  rates);  p  =  (cap(l), . . . ,  cap(n)),  where  cap(i)  is  the 
capacity2  of  the  ith  peer  (faster  peers  make  more  queries). 

•  (j  =  (<Ti, ,  crn),  the  relative  rates  at  which  the  peers  supply  items.  Reasonable 
value:  cr  =  r  (number  of  data  items  controlled  is  proportional  to  zone  size). 
Hot  spots  (Section  2. 2. 2. 3)  may  yield  a  cr  that  is  essentially  random  if  caching 
is  not  employed. 


We  will  employ  the  fact  that  the  weighted  geometric  mean 

VE?=i  «i 


wGM(r;  a)  =  I  ]^[ 


1=1 


(2.5) 


is  always  less  than  or  equal  to  the  weighted  arithmetic  mean 


wAM(r;  cl)  =  ^  (2.6) 

Z_/j=l  ai 

with  equality  iff  all  r\  are  equal.  The  proof  is  via  Jensen’s  inequality. 

It  will  be  useful  to  have  an  estimate  of  logfcwGM(r;  r ).  Let  f(x)  =  x\ogkx]  f 
is  infinitely  differentiable  and  concave  up  on  the  interval  /  =  [rmin,  rmax]  C  (0,  oo). 
Taylor’s  Theorem  says  that  for  any  points  a,  a  +  h  G  /  and  some  point  Gi  between  a 
and  a  +  h, 


f(a  +  h)  —  /(a)  +  f(a)h  +  ~f"{^h)h2 

X 

=  a  logfc  a  +  log k{ea)h  + 

>  a\ogka  +  \ogk(ea)h  + 

2  In  k  ^"max^  ' 

2A  peer’s  capacity  is  the  rate  at  which  it  can  process  P2P  traffic;  see  Section  2.3.1. 


(2.7) 

(2.8) 
(2.9) 
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The  mean  ^  Y^i=i  r*  =  so  let  s*  =  r*  —  1. 


Then  ]T"=1  =  0.  Finally, 


logfcWGM (r;  r) 


1 

-J2ri  l°^  ri 

1=1 

1  n 

-^/(i  +  Si) 

i=l 


n 


(o  +  logfc(e)s*  + 


1  _i  2\ 

21nA:rmaxS'V 


1 

2  In  k 
1 

2  In  k 


-i 


rmaxVar{'ri}- 


2. 2. 2. 2.1  Forward  query  path  lengths 

Theorem  2.2.  Consider  a  forward  de  Bruijn  DHT  with  n  peers. 

•  The  length  of  the  path  from  peer  i  to  any  given  destination  is 

if'd  =  [logfc  n  —  logfc  r f\ . 


•  Taking  a  weighted  average  over  all  peers  i  in  the  network  gives  the  the 
query  path  length 

Ifwd  =  logfc  n  -  logfcwGM(r;  p)  +  C 

for  some  0  <  C  <  1. 

Proof.  From  Section  2.2. 1.2.1,  the  length  of  the  path  from  i  to  any  b  is 

di  =  [log,,  (^i-1)]  =  logfc  n  -  logfc  n  +  ei} 
for  some  0  <  e%  <  1.  For  the  second  part  of  the  theorem,  we  have 

?wd  =  «(lo&  "  -  lo6t  r-  +  <0  J  f  Pi 

=  logt  n  -  logtwGM(r;  p)  +wAm(e;  {>) 


(2.10) 

(2.11) 

(2.12) 

(2.13) 

(2.14) 


(2.15) 
average 

(2.16) 

(2.17) 

(2.18) 
(2.19) 
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yielding  Eq.  (2.16). 


□ 


If  all  peers  have  the  same  request  rate  (p  =  1),  then 

logfcWGM (r;  1)  <  logfcwAM(r;  1)  =  log,,  1  =  0,  (2.20) 

so  the  average  query  path  length  is  at  least  logfc  n.  On  the  other  hand,  if  peers  with 
larger  zones  also  have  larger  p*  (e.g.,  p  =  (cap(l), . . . ,  cap(n ))  and  some  form  of  PLB 
has  been  performed),  then  wAM(r:  p)  >  1  and  the  average  query  path  length  can  be 
less  than  logfc  n. 


2. 2. 2. 2. 2  Reversed  query  path  lengths  We  begin  with  a  lemma. 


Lemma  2.1.  Consider  a  peer  j  in  a  reversed  de  Bruijn  DHT  with  n  peers.  If  b  G 
zone(j),  let  d b  be  minimal  such  that  ,bib2  ■  ■  ■  bdbc\c2  ■  ■  ■  G  zone(j)  for  all  c.  Then  the 
expected  value  of  db  over  zone(j)  is 


dj  =  logfc  n  -  logfc  rj  +  Cf  (2.21) 

where  0  <  C%- ]  <  2  +  ^  ^  <  4. 

Proof.  We  will  perform  the  computation  assuming  b  is  in  the  left  half  of  zoue(j); 
the  argument  for  the  right  half  is  symmetric.  Recall  from  Section  2.2. 1.2.2  that  db  is 
minimal  such  that  ( b  —  k~db,  b  +  k~db)  C  zone(j).  Thus  b  —  id(j )  G  [k~db,  k~db+l). 

Let  m  be  such  that  id(j)  +  [k~m,k~m+1)  contains  the  midpoint  of  zone(j);  we 
find  m  =  \—  logfc  y].  Then  db  is  a  step  function  on  its  half  of  zone(j),  given  by 


if  d  >  m  and  b  —  id(j )  G  [k  d,k  d+1), 


m  if  b  —  id(j)  G  [k  m,  ff  . 


(2.22) 
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The  expected  value  of  df,  is  found  by  weighted  average: 


dj  — 


J2  d(k-d+1~k~d) 


d=m-\- 1 


Zi 


+m{j~  k~m 

1,—m+l  -I  ~ 

rv  /On  .  /C" i 

- - 1-  m-1  — 

k- 1  2  '  2 


/- 
7  2 


=  [logfc  W  -  logfc  G  +  logfc  2]  + 


k 


— m+ 1 


(*-!)?■ 


(2.23) 

(2.24) 

(2.25) 

(2.26) 


We 


have  used  Gabriel’s  truncated  staircase,  X^m+i  ?'r*  =  (i_r)2  ^or  0  <  r  <  1- 


By  definition  of  m,  k  m+1  €  (y,/cyj.  Thus  the  final  term  satisfies 


fc  —  m+1 


(fcZi)  fezy]  G  (0,2],  yielding  Eq.  (2.21). 


□ 


We  now  proceed  to  the  theorem  that  is  the  main  result  of  this  subsection. 


Theorem  2.3.  Consider  a  reversed  de  Bruijn  DHT  with  n  peers. 

•  Averaging  over  all  destinations  b  chosen  uniformly  at  random  from  zone(j),  the 
expected  length  of  the  path  from  any  given  peer  to  b  is 

C”  ~  d  '  (logit  n  -  logfc  Tj  +  Cjf^  -  1,  (2.27) 

where  /i  <  4. 

•  Taking  a  weighted  average  over  all  destination  zones  j  in  the  network  gives  the 
the  average  query  path  length 


£rev  sa  p  ■  (logfc  n  -  logfcwGM(r;  cr)  +  Ck) .  (2.28) 


Here  C(f  ] 


and  C) t  are  positive  values  bounded  by  2  +  +  Wj-  <  4. 
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Proof.  In  Section  2.2. 1.2.2,  we  found  that  the  number  of  steps  £,1C'  =  2^]°B2  rE+l  —  1  < 
4 dj  —  1,  yielding  the  first  part  of  the  theorem  with  fi  =  4. 

For  the  second  part  of  the  theorem,  we  average  over  all  destination  zones  j, 

(n  \  n 

E  ff//ev  /  E  ai  <2-29> 

3= 1  /  3= 1 

n  \  —  1 

E  o  (2-3°) 

3=1  / 

=  p  ■  (logfc  n  —  logfcwGM(r;  a )  +wAm(C^;  ct))  (2-31) 

yielding  Eq.  (2.28).  □ 

In  practice,  n  will  depend  on  the  distribution  of  the  dj.  Letting  hj  =  [log2<ij], 

we  have  ii  =  2 hi+1/dj.  If  dj  were  uniformly  distributed  in  (2hj'_1,  2hj],  then  E(/i)  = 
4  In 2  «  2.8.  If  dj  were  distributed  as  l/x  in  (2hj~1,  2hj]  (such  as  by  Benford’s  Law), 
then  E(/i)  =  2/ In  2  «  2.9. 

If  received  requests  are  proportional  to  zone  size  (er  —  r),  then  logfcwGM(r;  r)  > 
0,  and  the  average  query  path  length  will  be  less  than  logfc  n.  In  particular,  we  saw 
above  that  logfcwGM(r;  r)  >  ypu; r'mIx^ar { r% } ?  so  high  variance  among  zone  sizes  has 
a  beneficial  effect  on  the  average  query  path  length. 

2. 2. 2. 3  Caching 

Dynamic  caching  is  an  important  tool  for  relieving  hot  spots  [28]  —  items  with  extreme 
popularity  that  cause  acute  congestion  for  the  peers  that  host  them.  Caching  is  a 
process  for  replicating  such  items  at  multiple  nodes,  to  distribute  the  workload. 


n 


E  '  (logfc  n  -  logfc  r 


+  C?h 


G=i 


Caching  is  compatible  with  PLB  and  provides  a  complementary  benefit.  Since 


PLB  is  based  on  the  actual  amount  of  traffic  processed  by  each  peer,  it  is  able  to  deal 
with  hot  spots  to  some  degree,  but  caching  can  be  faster  and  more  effective.  On  the 
other  hand,  caching  is  not  designed  to  alter  the  structure  of  the  network  in  order  to 
take  advantage  of  heterogeneity. 

2. 2. 2. 3.1  Forward  caching  Forward  de  Bruijn  DHTs  permit  a  simple  caching 
scheme  [17,  37].  A  hot  spot  at  location  b  e  [0,1)  should  be  replicated  at  the  k 
predecessor  locations  Pp  =  .xb\b2  •  •  • ,  x  G  7Lk.  Then  a  request  for  b  will  always 
pass  through  one  of  these  caches,  each  with  equal  probability.  The  replication  can  be 
repeated  as  necessary. 

2. 2. 2. 3. 2  Reversed  caching  The  analogous  scheme  does  not  work  for  reversed 
de  Bruijn  DHTs,  since  there  is  only  one  predecessor  location  kb  mod  1  to  the  hot 
spot  b.  A  routing-independent  caching  scheme  such  as  [28]  could  be  used.  Al¬ 
ternatively,  a  peer  that  routes  more  than  a  certain  number  of  requests  for  item  b 
can  begin  caching  the  item,  until  some  set  of  peers  together  controlling  the  inter¬ 
val  .6crit  •  •  ■  b2hbi  ■  ■  ■  b2h~i  ■  ■  ■  bib2bic ,  for  all  c  e  [0, 1),  distributes  the  workload  to  a 
satisfactory  level.  (Here  bcrit  is  the  first  value  for  which  the  workload  is  distributed 
satisfactorily.) 
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2.3  Proportionate  load  balancing 


In  this  section,  we  present  an  implementation  of  proportionate  load  balancing  in  the 
context  of  continuous-discrete  DHTs  (Section  1.3),  not  necessarily  de  Bruijn.  Using 
only  information  about  itself  and  the  two  peers  adjacent  to  it  in  the  keyspace  [0, 1), 
each  peer  will  increase  or  decrease  its  zone  size.  Assumption  1  indicates  that,  over 
time,  this  should  equalize  all  peers’  utilizations.  We  stress  that  in  simulations,  we 
will  have  all  peers  alter  their  zone  sizes  simultaneously  on  clock  ticks,  but  in  reality 
a  peer  can  alter  its  zone  size  more  or  less  often  depending  on  its  individual  needs. 

2.3.1  Definitions 

The  capacity  of  a  peer  (cap(i))  is  the  rate  at  which  it  is  able  to  process  traffic  for  the 
P2P  application.  Generally,  the  outbound  traffic  rate  will  be  the  limiting  factor  [29]. 
Plots  of  estimated  bottleneck  bandwidths  for  peers  in  the  Napster  and  Gnutella  [32] 
networks  are  given  in  Section  3.1  of  [46];  these  plots  show  that  there  are  frequently 
orders  of  magnitude  differences  between  peer  capacities. 

Capacity  is  difficult  to  measure  from  the  outside,  but  a  peer  can  keep  track  of 
its  own  capacity  with  relative  ease.  To  discourage  under-reporting  of  capacity,  peers 
might  limit  the  rate  at  which  they  transfer  hies  to  peer  i  to  be  some  multiple  of  i’s 
published  capacity.  This  does  not  affect  the  implementation  of  PLB.  Over-reporting  is 
a  more  complex  issue,  but  not  an  insurmountable  one.  The  predecessor  and  successor 
of  a  newly  joined  peer  can  restrict  the  growth  of  its  zone  by  modifying  the  algorithm 
of  Section  2. 3. 2. 2.  If  a  misbehaving  peer  does  obtain  a  large  zone  and  fails  to  pass  on 
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all  its  messages,  it  can  be  treated  like  a  failed  peer  as  in  Section  2.3.2. 1. 

The  load  of  a  peer  (load(i))  is  the  rate  at  which  the  network  asks  it  to  pro¬ 
cess  P2P  traffic.  A  peer  can  keep  track  of  this  value,  perhaps  as  an  exponentially 
weighted  moving  average.  Load  is  a  dynamic  quantity  that  can  depend  on  relative 
zone  size  (Theorem  2.1),  global  network  structure  (Fig.  2.3),  and  key  popularity  (Sec¬ 
tion  2. 2. 2. 3)  in  a  complicated  way. 

The  utilization  of  a  peer  ( util(i ))  is  a  dimensionless  quantity  equal  to  the  ratio 
of  its  load  and  capacity:  util(i)  =  In  order  for  a  peer  to  keep  up  with  its 

requests,  its  utilization  should  remain  below  1. 

2.3.2  Implementation 

In  our  implementation,  peer  i  is  only  required  to  keep  track  of  limited  information: 

For  the  DHT:  Its  zone  zone(i )  =  [id(i),  id(succ(i)))  and  the  data  items  that 
map  to  it.  Links  to  its  neighbors  in  the  underlying  graph,  for  routing.  Links  to  peers 
pred(i)  and  succ(i),  for  consistency. 

For  PLB:  The  capacity,  load,  and  utilization  of  itself  and  peers  pred(i )  and 
succ(i). 

2. 3. 2.1  Arriving  and  departing  peers  (churn) 

When  a  new  peer  wants  to  join  a  P2P  network,  typically  it  must  know  how  to  contact 
some  existing  peer,  who  helps  the  new  peer  find  its  proper  neighbors.  Several  methods 
for  assigning  a  zone  to  a  new  peer  were  analyzed  in  [54],  In  single-point  random- split, 
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a  randomly  chosen  point  in  the  keyspace  becomes  the  new  peer’s  id.  In  single-point 
center-split,  the  peer  that  controls  a  randomly  chosen  point  gives  exactly  half  of  its 
zone  to  the  new  peer.  In  multi-point  center-split  (random  or  semi-deterministic),  a 
fixed  number  of  points  are  sampled;  the  one  belonging  to  the  largest  zone  becomes 
the  target  of  a  center-split.  As  an  example  of  the  analyses,  we  have  that  under 
single-point  center-split  the  largest  and  smallest  zones  are  fl(logn)  times  larger  and 
Q  (V21og2")  times  smaller  than  average,  with  high  probability. 

The  center-split  methods  are  readily  generalizable  to  proportional- split  methods, 
where  the  target  zone  is  split,  not  in  half,  but  in  proportion  to  the  capacities  of  the 
new  and  existing  peers3.  The  intent  is  to  equalize  utilizations  rather  than  zone  sizes. 
In  a  multi-point  method,  the  point  belonging  to  the  peer  with  the  highest  utilization 
becomes  the  target. 

In  our  simulations,  we  will  use  single-point  proportional-split,  which  works  as 
follows:  If  the  network  has  n  —  1  peers,  then  the  nth  peer  chooses  a  random  point 
id  G  1C,  which  will  lie  in  the  zone  of  some  existing  peer  i.  Peer  n  then  receives  the 
identifier  id(n )  G  zone(i)  such  that  f*-  =  (Alternatively,  id(n )  could  be  chosen 

based  on  the  load  distribution  function  gi{x)  discussed  in  Section  2. 3. 2. 2. 2.) 

When  peer  i  leaves  the  P2P  network,  its  adjacent  peers  pred(i )  and  succ(i) 

should  take  control  of  its  zone.  If  peer  i  announces  its  departure,  the  adjacent  peers 

can  split  zone{i )  in  proportion  to  their  capacities,  and  take  responsibility  for  the  data 

items  previously  managed  by  i.  If  peer  i  fails  or  simply  drops  out,  some  information 
3More  generally,  it  may  be  split  in  proportion  to  some  function  of  the  capacities  of  the  new  and 
existing  peers,  depending  on  the  nature  of  the  relationship  in  Assumption  1. 
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may  be  temporarily  lost.  The  Chord  coping  mechanism  [49]  is  applicable  to  our 
DHTs;  an  analysis  of  coping  mechanisms  is  beyond  the  scope  of  this  paper.  Typically, 
periodic  re-publication  of  data  items  is  required. 

2. 3. 2. 2  Existing  peers 

Proportional-split  joining  methods  are  insufficient  to  achieve  balanced  utilizations. 
Part  of  this  is  due  to  the  unpredictable  nature  of  peer  departures,  but  at  least  as 
important  is  the  fact  that  local  DHT  traffic  load  depends  on  the  global  structure  of 
the  network.  Fig.  2.3  shows  the  distribution  of  load  across  a  particular  peer’s  zone 
for  an  actual  simulation,  a  reversed  de  Bruijn  network  in  which  the  request  rates  p 
(Section  2. 2. 2. 2)  were  all  equal  and  destinations  b  E  1C  were  chosen  with  uniform 
probability. 

The  dependence  of  load  on  network  structure  means  that  it  is  not  practical 
(or  perhaps  even  possible)  to  find  an  exact  solution  to  the  problem  of  utilization 
balancing,  even  if  the  set  of  peers  is  held  constant.  This  is  in  contrast  to  processor 
networks  or  parallel  databases,  where  an  exact  solution  can  be  found  as  the  result  of  a 
global  linear  equation,  and  the  challenge  is  to  find  locally-computable  approximations 
to  this  global  solution. 

Our  implementation  of  PLB  has  similarities  to  the  NbrAdjust  operation  of 
[18]  and  the  first-order  diffusion  schemes  described  in  [14],  but  we  also  address  the 
non-uniformity  depicted  in  Fig.  2.3,  and  we  incorporate  fail-safes  to  make  sure  that 
zone  sizes  do  not  change  too  quickly.  Pseudocode  for  the  local  implementation  of 
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g.(x):  Fraction  of  peer  i’s  traffic 
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Figure  2.3:  Traffic  load  is  not  distributed  uniformly  within  zones,  even  under 
uniform  selection  of  sources  and  destinations.  This  graph  depicts  the  distri¬ 
bution  of  load  across  a  particular  peer’s  zone  in  a  reversed  de  Bruijn  network. 
In  the  simulation,  19%  of  the  traffic  went  through  the  second  portion  of  the 
zone,  and  only  6%  through  the  fifth  portion.  Call  this  step  function  gi{x)  for 
peer  i. 
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PLB  is  given  in  Fig.  2.4. 


2. 3. 2. 2.1  The  Balance  operation  Lines  1  and  2  of  Balance  assign  h  to  be 
the  peer  with  higher  utilization,  and  £  the  peer  with  lower  utilization.  Line  3  tests 
whether  util(h )  is  sufficiently  higher  than  util(£)  for  the  GiveZone  operation  to  be 
performed  in  line  4. 

PLB  requires  two  tolerance  values,  tm in,  td iff  £  [0, 1).  Their  meaning  is  as  follows. 
If  util(£ )  =  0,  we  change  id(i)  if  util(h)  —  util{£)  >  fmin.  On  the  other  hand,  if  the 
util(£ )  =  1,  we  change  id(i )  if  util(h )  —  util(£)  >  fmin  +  Giff-  For  other  values  of 
util(£),  we  interpolate  linearly.  Smaller  values  of  tmin  and  td iff  cause  the  network  to 
take  longer  to  settle  down,  but  a  more  even  distribution  of  utilizations  is  achieved. 

2. 3. 2. 2. 2  The  GiveZone  operation  Line  1  of  GiveZone  determines  the  frac¬ 
tion  of  peer  K s  load  that  should  be  transferred  to  peer  £  so  that  util(h)'  =  util(£)'. 
Solving  the  equation  for  /  yields  /  =  ^ g^>.  Be- 
cause  i  communicates  with  both  pred(i)  and  succ(i),  zone(i)  may  be  expanded  or  con¬ 
tracted  on  both  ends  simultaneously,  so  we  use  the  modified  value  loadfrac(h,  £)  =  \f. 

Recall  that  the  distribution  of  traffic  load  is  nonuniform  within  zones.  In  line 
2,  peer  h  determines  the  fraction  zonefrac(h,  £)  of  its  zone  (adjacent  to  id(i))  that 
corresponds  to  a  ■  loadfrac(h,  £)  of  its  load.  (The  parameter  a  G  (0, 1]  is  tunable;  we 
use  a  =  1.)  In  our  simulations,  we  had  each  peer  keep  a  histogram  of  the  activity  in 
s  =  10  equal  sub-intervals  of  its  zone,  generating  a  step  function  fjh  like  in  Fig.  2.3. 
If  h  =  i,  the  integral  is  Gh,e(x)  =  gh{t)dt.  If  h  =  pred(i),  the  integral  is  Ghj(x)  = 
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Balance  (peer  i) 

1  Let  h  <—  argmax{ util(i),  util(pred(i ))} 

2  Let  £  argminj util(i),  util(pred(i ))} 

3  IF  util(h)  >  (1  +  tdis)  util(£)  +  tmin  THEN 

4  GiveZone(/?,,  £) 


GiveZone  (peer  h,  peer  £) 

1  Let  loadfrac(h,  £)^\- 

2  Let  zonefrac(h,  £)  G^\(aloadfrac(h,  £)) 

3  Let  r  <—  min  ^zonefrac(h,  £),  | , 

4  Adjust  id(i)  by  r  •  ^ 


Figure  2.4:  Pseudocode  for  PLB  (Section  2. 3. 2. 2).  Balance  describes  the 
conditions  under  which  PLB  is  performed.  Here  fmin,  Liiff  £  [0, 1).  GiveZone 
determines  the  new  boundary  id(i)/  between  zone(i)  and  zone(pred(i)) .  In 
simulations,  we  have  all  peers  execute  Balance  simultaneously  on  clock  ticks, 
but  in  reality  it  can  be  executed  at  any  time  depending  on  each  peer’s  needs. 
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x  =  fraction  of  zone(i) 


id(i-l) 


id(i)’ 


id(i+1) 


(a)  Contracting  zone(i) 


x  =  fraction  of  zone(i-l) 


(b)  Expanding  zone(i) 


Figure  2.5:  In  Fig.  2.5(a),  peer  i  contracts  its  zone  by  increasing  id(i)  to  id(i)' . 

In  Fig.  2.5(b),  peer  i  expands  its  zone  by  decreasing  id(i)  to  id(i)'. 

Ji_  Qh(t)dt.  Notice  that  is  easy  to  calculate  since  gu  is  a  step  function. 

Line  3  ensures  that  z'h  =  Zh  —  TZh  >  \zh  and  z'f  —  +  TZh  <  2 Z£,  so  the  zone 

sizes  do  not  adjust  too  quickly.  It  is  easy  to  repeat  the  Balance  operation,  but  it  is 
difficult  to  recover  from  a  suddenly  over-utilized  peer. 

In  line  4,  id(i )  is  adjusted  by  the  proper  fraction  of  Zh-  It  is  increased  if  h  =  i, 
and  decreased  if  h  =  pred(i).  See  Fig.  2.5. 
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2.3.3  Data  and  rewiring  costs 


When  peer  i  changes  its  id(i)  based  on  PLB,  the  interval  between  the  old  id(i )  and 
the  new  id(i)'  changes  ownership  as  discussed  in  the  last  section.  The  previous  owner 
is  h  (i  or  pred(i),  the  peer  with  higher  utilization),  while  the  new  owner  is  the  other 
peer  £.  There  is  a  data  cost  associated  with  moving  information  about  data  items 
from  h  to  £  to  maintain  the  proper  definition  of  their  zones.  This  cost  is  proportional 
to  the  relative  size  of  the  interval,  rr^,  that  was  transferred. 

There  is  also  a  rewiring  cost  for  adjusting  neighbor  connections  maintain  the 
structure  of  the  DHT.  When  id(i )  is  changed,  links  to/from  peer  h  may  be  removed, 
and  links  to/from  peer  £  may  be  added.  For  each  link  involving  h,  we  check  whether 
the  link  should  be  removed  from  h,  and  whether  a  corresponding  link  should  be  added 
for  £.  (These  are  separate  questions  because  both  h  and  £  may  simultaneously  link  to 
the  same  peer.)  Any  links  not  involving  h  will  be  unchanged.  A  special  case  arises  if 
h  has  a  (virtual)  link  to  itself;  then  we  must  check  all  four  combinations  of  old  and 
new  interval  owners:  h  — >  h,  £  — >  h,  h  — >  £,  and  £  — ■»  £. 

We  summarize  these  results  in  a  theorem. 

Theorem  2.4.  Consider  continuous- discrete  DHTs  as  in  Section  1.3.  Suppose  id(i ) 
is  adjusted  by  t  ■  Zh  in  the  GiveZone  operation. 

•  Averaging  over  all  distributions  of  data  items,  the  expected  data  cost  is 

TZh  D^j  =  rrhD.  (2.32) 

Here  Dj  is  the  number  of  data  items  stored  by  peer  j ,  and  D  —  ^  Y^j=i  Dj  is 
the  average  number  of  stored  (and  thus  of  shared)  data  items  per  peer. 
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•  The  rewiring  cost  is  at  most  2(out-degree  of  h )  +  2 {in-degree  of  h )  +  4. 


39 


2.4  Simulation  experiments 


We  now  use  simulations  to  evaluate  the  performance  of  our  proportionate  load  bal¬ 
ancing  implementation.  These  simulations  are  all  of  P2P  networks  based  on  various 
types  of  de  Bruijn  DHTs. 

One  key  feature  of  a  network,  which  we  will  report  on  for  each  experiment,  is  its 
maximum  utilization.  The  maximum  utilization  is  the  greatest  utilization  experienced 
by  any  peer  in  the  network.  A  large  maximum  utilization  (near  1)  implies  that  some 
peer  is  near  capacity,  and  will  start  dropping  queries  if  the  request  rates  increase.  A 
small  maximum  utilization  (near  0)  implies  that  all  peers  are  experiencing  relatively 
low  load,  and  the  request  rate  can  be  increased  substantially  without  overloading  the 
network. 

2.4.1  Setup 

We  perform  iterations  of  PLB  on  static  and  dynamic  networks  of  various  sizes,  ranging 
from  64  to  8192  peers.  During  one  iteration ,  all  utilizations  are  computed,  and  then 
each  peer  performs  the  Balance  operation  (Section  2. 3. 2. 2)  once. 

The  capacity  of  each  peer  is  assumed  to  be  constant  during  the  simulation.  Ca¬ 
pacities  are  selected  from  a  probability  distribution  that  approximates  the  estimated 
bottleneck  upstream  bandwidths  of  peers  in  the  Gnutella  network  [46].  This  distri¬ 
bution  is  described  in  Table  2.1.  Capacities  selected  from  this  distribution  will  be 
highly  variable,  indicated  by  the  standard  deviation  of  F  being  larger  than  the  mean. 

The  request  rate  pt  of  each  peer  (Section  2. 2. 2. 2)  is  also  assumed  to  be  constant 
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Table  2.1:  Statistical  properties  of  our  capacity  distribution. 


Property 

Value 

CDF 

F(x)  =  2^3  log10  § 

Mean 

940  Kbps 

Standard  deviation 

1.5  Mbps 

Median 

350  Kbps 

Range 

25  Kbps  <  x  <  5  Mbps 

during  the  simulation.  We  test  two  reasonable  values  of  p:  1  and  (cap(l), . . . ,  cap(n)) 
(Section  2. 2. 2. 2).  In  the  simulations  presented  here,  the  supply  rate  a,  of  each  peer 
is  taken  to  be  proportional  to  its  relative  zone  size  r*.4 

We  estimate  the  utilization  of  each  peer  by  tracing  paths  from  random  source 
peers  to  random  destination  keys.  The  number  of  times  a  peer  is  traversed  during 
this  process  allows  us  to  estimate  its  relative  traffic  load.  An  average  of  100  paths  per 
source  peer  were  traced;  results  were  similar  for  different  choices  of  the  traced  paths. 

We  tested  various  protocol  parameters,  and  found  no  unpredictable  results.  The 
parameters  used  herein  are  as  follows,  de  Bruijn  graphs:  k  =  8,  D  =  8.  BALANCE 
operation:  tmin  =  0.03,  tdiff  =  0.07. 

2.4.2  Static  networks 

In  our  static  simulations,  no  peers  joined  or  left  the  system.  We  assigned  capacities 

to  the  n  peers  and  fixed  their  order  in  the  keyspace.  We  report  on  four  stages  of  PLB. 

In  each  we  measured  network  features  such  as  maximum  utilization,  mean  utilization, 
Simulations  with  smooth  but  very  nonuniform  distributions  of  data  items  showed  that  PLB  copes 
well  with  different  values  of  er;  space  prohibits  their  presentation  here.  Caching  can  be  employed  to 
level  isolated  hot  spots,  as  described  in  Section  2. 2. 2. 3. 
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and  average  query  path  length. 


0.  As  a  baseline,  we  tested  the  uniform  network  (all  peers  have  Zi  —  -). 

1.  We  then  rebuilt  the  network  with  the  single-point  proportional-split  method  of 
Section  2.3.2. 1. 

2.  We  performed  five  iterations  of  PLB  on  this  network  (to  evaluate  realistic  per¬ 
formance). 

3.  We  continued  performing  500  iterations  for  some  simulations  (to  evaluate  lim¬ 
iting  behavior).  The  worst-case  value  of  each  feature  between  450  and  500 
iterations  was  used. 

2. 4. 2.1  Uniform  request  rate 

In  these  simulations,  p  oc  1.  At  each  network  size,  30  simulations  were  performed  for 
stages  1  and  2.  Because  of  time  considerations,  only  6  simulations  were  continued  to 
stage  3. 

Figs.  2.6  and  2.7  show  how  the  maximum  utilization  drops  under  PLB.  The 
inverse  of  the  maximum  utilization  determines  how  much  traffic  the  network  can 
handle  without  any  peer  becoming  over-utilized. 

Relative  to  the  corresponding  uniform  networks,  forward  de  Bruijn  DHTs  saw 
a  greater  decrease  in  maximum  utilization  in  stages  1  and  2  than  did  reversed  de 
Bruijn  DHTs.  However,  the  starting  values  in  stage  0  were  about  twice  as  great.  The 
limiting  behavior  of  both  cases  was  about  the  same,  with  the  maximum  utilization 
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Ratio  to  max.  util, 
in  forward  uniform  network 


Figure  2.6:  Results  of  the  static  simulations  in  Section  2. 4. 2.1,  with  uniform 
request  rate  p  =  1.  The  results  for  forward  de  Bruijn  networks  are  depicted 
here.  Multiple  simulations  were  performed.  These  plots  show  median  values 
connected  by  lines,  with  boxes  delimiting  the  25th  and  75th  percentiles,  and 
whiskers  showing  the  full  extent  of  the  data.  Figs.  2.6  and  2.7  show  how 
the  maximum  peer  utilization  decreases  under  PLB.  All  measurements  are 
relative  to  the  maximum  utilization  in  the  corresponding  uniform  network 
(peer  capacities  and  ordering  are  the  same  as  the  original  network,  but  all 
zones  are  given  equal  size).  The  solid  line  shows  networks  built  with  single¬ 
point  proportional-split  (Section  2. 3. 2.1).  The  dashed  line  shows  the  same 
proportional-split  networks  after  5  iterations  of  Balance.  The  dotted  line 
shows  the  limiting  behavior  of  PLB,  after  450  iterations. 
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Ratio  to  max.  util, 
in  reverse  uniform  network 
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Figure  2.7:  Results  of  the  static  simulations  in  Section  2. 4. 2.1,  for  reversed  de 
Bruijn  networks.  This  figure  complements  Figure  2.6,  which  shows  the  results 
for  forward  de  Bruijn  networks. 
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Figure  2.8:  Sample  result  of  the  static  simulations  in  Section  2. 4. 2.1.  For  a 
fixed  injection  rate  p,  we  compare  the  maximum  utilization  for  forward  and 
reversed  de  Bruijn  DHTs  as  PLB  is  performed.  The  graph  depicts  the  results 
of  one  simulation  of  each,  with  2048  peers.  Local  fluctuations  are  much  higher 
in  the  forward  case. 

being  5%  to  10%  of  its  original  value. 

As  predicted  in  Section  2. 2. 2. 2,  query  path  lengths  decreased  (by  10%)  for 
reversed  de  Bruijn  DHTs,  and  increased  (by  10%)  for  forward  de  Bruijn  DHTs.  In 
stage  0,  the  reversed  query  path  lengths  were  about  70%  longer  than  the  forward 
query  path  lengths. 
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2. 4. 2. 1.1  Volatility  Figure  2.8  compares  the  actual  values  of  the  maximum  uti¬ 
lization  for  forward  and  reversed  de  Bruijn  DHTs.  In  the  limit  of  a  large  number  of 
iterations,  their  best-case  behavior  is  similar,  but  oscillations  for  the  forward  DHTs 
yield  a  higher  worst-case  behavior. 

The  sudden  spikes  in  utilization  seen  in  Figure  2.8  most  often  occur  when  a 
low-capacity  peer  becomes  overwhelmed,  usually  for  one  of  two  reasons,  (i)  The  peer 
might  expand  its  zone  and  acquire  more  traffic  than  expected,  due  to  imbalances  in  the 
keyspace.  The  imbalances  are  more  profound  in  forward  de  Bruijn  graphs,  because  of 
their  more  primitive  routing  procedure,  (ii)  A  different  peer  might  change  the  size  of 
its  zone,  thus  forcing  it  to  connect  to  the  peer  in  question.  This  is  more  problematic  in 
the  forward  case  because  a  small  peer  is  likely  to  have  only  one  outgoing  connection, 
thus  overwhelming  its  sole  neighbor. 

2. 4. 2. 2  Varied  request  rate 

In  these  simulations,  p  oc  (cap(l), . . . ,  cap(n)).  At  each  network  size,  30  simulations 
were  performed  for  stages  1  and  2.  Of  these,  6  were  continued  to  stage  3.  See 
Figure  2.9. 

Results  for  maximum  utilization  were  very  similar  to  the  simultations  with 
uniform  rerquest  rate,  as  seen  in  Figure  2.9.  Forward  de  Bruijn  DHTs  again  saw  a 
greater  decrease  in  maximum  utilization  in  stages  1  and  2  than  did  reversed  de  Bruijn 
DHTs.  The  starting  values  in  stage  0  were  again  about  twice  as  great.  However,  the 
limiting  behavior  in  forward  DHTs  was  slightly  better  than  reversed. 
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Number  of  peers  (log  scale) 


(a)  Forward  de  Bruijn 


(b)  Reversed  de  Bruijn 


Figure  2.9:  Results  of  the  static  simulations  in  Section  2. 4. 2. 2,  with  capacity- 


proportional  request  rate.  See  Figure  2.6  for  an  explanation. 
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Query  path  lengths  again  decreased  for  reversed  de  Bruijn  DHTs.  Forward  de 
Bruijn  DHTs  saw  a  very  slight  decrease.  In  stage  0,  the  reversed  query  path  lengths 
were  about  65%  longer  than  the  forward  query  path  lengths. 

2.4.3  Dynamic  networks 

These  simulations  model  networks  that  have  reached  a  steady-state  size,  but  with  a 
fluctuating  (dynamic)  population  of  peers.  We  start  with  the  average  network  size 
n  =  512,  and  allow  peers  to  join  and  leave  the  system  according  to  the  rules  of 
Section  2.3.2. 1. 

Fix  a  churn  rate  y.  Between  each  iteration  of  PLB,  we  independently  allow 
each  peer  to  drop  out  of  the  network  with  probability  y,  and  we  also  insert  512y  new 
peers  with  capacities  chosen  from  the  distribution  in  Table  2.1.  This  process  will  tend 
to  keep  the  network  size  close  to  512. 5  How  fast  do  peers  come  and  go  in  reality? 
A  detailed  study  of  churn  rates  ([50],  Figure  6)  reveals  that  in  the  real-world  DHT 
Kad,  approximately  80%  of  the  peers  in  the  network  at  any  given  time  have  been 
present  for  at  least  one  hour,  while  approximately  95%  have  been  present  for  at  least 
ten  minutes.  Therefore,  only  a  small  percentage  of  the  peers  will  have  changed  if  we 
execute  PLB  every  few  minutes. 

We  executed  12  simulations  of  500  iterations  for  forward  and  reversed  networks 
for  values  of  y  =  2“15,  2~14, . . . ,  2~3.  The  results  for  small  y  were  similar  to  Figure  2.8, 

as  expected.  A  representative  example  for  larger  y  appears  in  Figure  2.10. 

5If  the  network  size  is  rij  after  i  iterations,  then  E  (| rij+i  —  512|  |  rij)  =  |(rij  —  rux  +  512y)  —  512|  = 
(1  —  X)\ni  —  512 1  <  |ru-512|. 
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Util,  relative  to  forwarduniform  network 


Figure  2.10:  Sample  result  of  the  dynamic  simulations  in  Section  2.4.3  (for¬ 
ward  network  pictured).  The  churn  rate  is  x  =  2~5  ~  3.1%.  For  clarity,  only 
200  iterations  are  displayed.  Forward  and  reversed  networks  displayed  statis¬ 
tically  similar  reactions  to  churn  rate.  We  see  that  the  maximum  utilization 
is  highly  volatile,  but  the  99th  percentile  of  utilizations  is  relatively  stable, 
indicating  that  the  vast  majority  of  peers  in  dynamic  networks  will  benefit 
from  PLB. 
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For  definiteness,  we  define  a  spike  as  an  increase  in  one  step  by  at  least  a  factor 
of  two,  and  the  spike  recovery  time  as  the  number  of  steps  to  return  to  within  10% 
of  the  pre-spike  value.  First  consider  the  maximum  peer  utilization,  as  we  have  in 
other  experiments.  As  x  becomes  large  (we  tested  up  to  x  —  2-3  ~  13%),  the 
median  distance  between  spikes  and  the  median  spike  recovery  time  both  approach 
approximately  5  steps.  Now  consider  the  99th  percentile  of  the  peer  utilizations, 
depicted  by  the  bottom  line  in  Figure  2.10.  Ignoring  a  very  few  highly-affected  peers 
increases  the  median  distance  between  spikes  to  approximately  75,  without  greatly 
affecting  the  spike  recovery  time. 

Note  that  we  are  assuming  no  correlation  between  session  length  and  capacity. 
It  is  plausible  to  think  that  in  reality,  peers  with  short  session  times  also  have  low 
capacities.  Such  behavior  would  reduce  the  fluctuations  in  maximum  utilization,  since 
the  removal  of  a  peer  with  a  small  zone  has  less  impact  than  the  removal  of  one  with 
a  large  zone. 

2.4.4  Hierarchical  DHTs 

Hierarchical  P2P  systems  [19]  divide  their  peers  into  two  (or  sometimes  more)  classes. 
One  class  consists  of  superpeers ,  with  high  capacity  and  stability,  and  one  consists 
of  regular,  transient  peers.  We  implemented  a  simple  hierarchical  system  based  on 
reversed  de  Bruijn  graphs,  as  follows. 

Given  a  threshold  capacity  c,  we  label  peer  i  a  superpeer  if  cap(i)  >  c,  or  a 
regular  peer  otherwise.  We  form  a  de  Bruijn  DHT  from  the  superpeers;  then  each  of 
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Ratio  to  max.  util, 
in  forward  uniform  flat  network 
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Figure  2.11:  Results  of  the  hierarchical  simulations  in  Section  2.4.4.  The 
results  for  forward  de  Bruijn  neigworks  are  depicted  here;  see  Figure  2.12 
for  the  results  for  reversed  de  Bruijn  networks.  Peers  were  grouped  by  the 
threshold  capacity  cmean  =  940  Kbps,  and  results  are  shown  relative  to  the 
corresponding  flat  uniform  network  (no  hierarchy,  all  zones  equal  size).  The 
solid  line  depicts  a  hierarchical  uniform  graph,  while  the  dashed  line  shows  a 
hierarchical  graph  to  which  PLB  has  been  applied  for  5  iterations. 
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Ratio  to  max.  util, 
in  reverse  uniform  flat  network 
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Figure  2.12:  Results  of  the  hierarchical  simulations  in  Section  2.4.4,  for  re¬ 
versed  de  Bruijn  networks.  This  figure  complements  Figure  2.11,  which  shows 
the  results  for  forward  de  Bruijn  networks 
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the  regular  peers  is  assigned  to  a  superpeer  at  random.  All  peers  share  data,  but  a 
superpeer  handles  all  of  the  routing  of  requests  for  its  associated  regular  peers. 

When  building  the  network  with  the  single-point  proportional-split  method 
(Section  2.3.2. 1),  we  split  zones  in  proportion  to  the  square  roots  of  the  capacities 
of  the  new  and  existing  superpeers.  This  is  because  the  expected  amount  of  traffic 
through  a  superpeer  depends  both  on  its  zone  size,  and  on  the  number  of  regular 
peers  attached  to  it,  which  is  proportional  to  its  zone  size.  The  BALANCE  operation 
was  unchanged. 

Figures  2.11  and  2.12  show  results  for  uniform  hierarchical  networks  and  PLB- 
iterated  hierarchical  networks,  both  compared  to  uniform  non-hierarchical  (flat)  net¬ 
works.  We  note  two  features  of  these  plots.  First,  uniform  hierarchical  networks 
have  maximum  utilizations  that  are  about  3%  (forward)  to  6%  (reversed)  of  the  value 
for  a  flat  network,  which  is  a  large  improvement.  Second,  executing  PLB  on  these 
hierarchical  networks  yields  maximum  utilizations  that  are  as  good  as  or  better  than 
those  for  the  uniform  hierarchical  networks. 
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2.5  Discussion 


We  have  presented  the  mechanics  of  forward  and  reversed  de  Bruijn  DHTs,  under 
the  continuous-discrete  framework.  We  analyzed  the  effect  of  relative  zone  size  on 
expected  degree  in  each  type  of  DHT.  Furthermore,  we  investigated  the  effect  of  zone 
size  distribution  on  the  average  query  path  length.  Theorems  2.1,  2.2,  and  2.3  of 
Section  2.2.2  describe  the  beneficial  effect  heterogeneity  can  have. 

We  have  also  presented  a  proportionate  load  balancing  algorithm,  which  it¬ 
eratively  optimizes  continuous-discrete  DHTs  by  equalizing  peer  utilizations.  After 
analyzing  the  data  and  rewiring  costs  involved,  we  performed  simulations  of  static, 
dynamic,  and  hierarchical  networks  to  verify  the  efficacy  of  PLB  in  reducing  strain 
on  peers.  We  found  that  the  maximum  relative  load  on  any  peer  dropped  to  approx¬ 
imately  20%  of  the  value  for  a  uniformly-structured  network  after  just  five  iterations 
of  PLB,  and  to  less  than  10%  in  the  limit  of  450-500  iterations.  A  more  complete 
practical  algorithm  could  be  a  subject  of  future  research. 

Networks  built  with  PLB  offer  a  continuum  of  peer  roles,  from  server-like  to 
client-like,  so  each  peer  can  choose  an  appropriate  role  for  itself  depending  on  the 
current  makeup  of  the  network.  This  differs  from  hierarchical  schemes  [19]  which 
treat  some  peers  almost  exclusively  as  servers  and  the  others  only  as  clients. 

We  found  that  forward  de  Bruijn  DHTs  are  more  volatile  in  their  response  to 
PLB  than  reversed  de  Bruijn  DHTs  (Section  2. 4. 2. 1.1).  They  also  experienced  larger 
starting  values  of  maximum  utilization,  by  approximately  a  factor  of  two.  Although 
forward  DHTs  have  shorter  average  query  path  length  in  uniformly-structured  net- 
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works,  under  heterogeneity  we  have  that  the  query  path  length  may  increase  for 
forward  DHTs,  while  for  reversed  DHTs  it  decreases  by  rmL Var { r, } .  This  is 

because  in  the  forward  case,  a  large  zone  size  only  benefits  the  owner,  while  in  the  re¬ 
versed  case,  a  large  zone  size  benefits  all  the  peers.  We  believe  these  features  present 
a  strong  case  for  the  consideration  of  reversed  de  Bruijn  DHTs. 
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Chapter  3 


Identification  of  peer-to-peer  hosts 
3.1  Introduction 

Universities  and  corporations  have  a  particular  interest  in  identifying  peer-to-peer 
traffic  in  their  networks.  Their  concerns  are  twofold:  bandwidth  and  security.  With 
regard  to  bandwidth,  peers  in  peer-to-peer  networks  typically  upload  more  than 
clients  in  client-server  systems  [45],  which  puts  strain  on  the  networks  they  inhabit. 
With  regard  to  security,  sharing  pirated  content  can  result  in  costly  legal  action,  or 
applications  and  content  may  contain  viruses,  spyware,  or  other  malware  [7]. 

Our  focus  on  this  chapter  is  on  security  concerns,  so  our  goal  is  to  identify  the 
offending/vulnerable  hosts.  In  contrast,  when  the  major  concern  is  for  bandwidth, 
the  goal  is  to  identify  individual  traffic  flows,  or  even  packets.  Higher  accuracy  can 
be  achieved  when  identifying  hosts,  because  more  information  is  available:  we  can 
consider  all  the  flows  a  host  is  involved  in.  Additionally,  network  monitoring  is  work¬ 
intensive.  We  are  able  to  achieve  acceptable  results  with  NetFlow-stylc  [8]  data, 
while  effective  flow-identification  schemes  tend  to  require  packet-level  information 
[3,  55,  35], 

We  will  use  behavioral  commonalities  to  identify  peer-to-peer  applications.  Our 
data  will  be  IP  traffic  passively  observed  at  set  monitoring  locations.  Examples  of 
observable  behaviors  are  average  flow  size,  or  number  of  IP  addresses  communicated 
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with  in  a  certain  period  of  time.  Message  passing,  measured  in  “triples,”  will  be  a 
particularly  important  behavior.  We  will  intentionally  ignore  other  behaviors,  such  as 
port  information  [34],  or  the  specific  timing  of  requests,  because  they  are  application- 
specific.  We  will  also  attempt  to  make  the  identification  as  close  to  real-time  as 
possible,  rather  than  performing  extensive  post-processing  analysis  of  network  struc¬ 
ture. 

3.1.1  Related  work 

Most  of  the  literature  on  identifying  peer-to-peer  traffic  has  focused  on  labeling  flows 
(or  IP  pairs)  rather  than  individual  hosts.  When  applicable,  payload  examination  can 
be  a  very  effective  method  of  classifying  particular  P2P  applications  [47,  1].  However, 
encryption  and  data  volume  are  increasingly  rendering  payload  methods  inapplicable. 

It  is  also  possible  to  identify  flows  based  only  on  features  extracted  from  meta¬ 
data.  These  features  are  typically  processed  by  machine  learning  classification  tech¬ 
niques.  Some  studies  have  attempted  to  identify  specific  applications,  such  as  Kazaa 
or  eDonkey.  Two  such  studies  achieved  80-90%  accuracy,  one  [3]  using  only  the  sizes 
of  its  first  several  data  packets  in  the  flow,  and  another  [55]  using  a  wider  variety  of 
features,  including  packet  interarrival  times. 

Other  studies,  which  are  more  in  line  with  our  goals,  have  attempted  to  identify 
the  categories  of  traffic  that  a  flow  belongs  to,  such  as  P2P  or  email.  One  such 
study  [35]  used  features  including  packet  interarrival  times  and  TCP  ports  to  achieve 
up  to  95%  accuracy  for  certain  application  categories,  though  only  up  to  55%  for 
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P2P.  Another  [26]  used  features  including  port  counts  to  identify  99%  of  P2P  flows 
with  a  10%  false  positive  rate,  but  used  a  fixed  algorithm  rather  than  adaptable 
machine  learning  techniques.  Though  these  studies  did  not  rely  on  prespecified  port 
information,  the  use  of  port  numbers  at  all  is  increasingly  unreliable,  because  it  is 
possible  for  users  to  employ  random  port  numbers. 

A  smaller  number  of  studies  have  discussed  identification  of  individual  hosts, 
rather  than  flows.  BLINC  [27]  used  features  including  port  counts  to  classify  hosts  into 
application  categories,  then  used  that  information  to  classify  the  corresponding  flows, 
achieving  95%  accuracy  on  80-90%  of  the  flows.  Some  others  focused  on  more  narrow 
subsets  of  P2P  applications:  one  [48]  primarily  investigated  eDonkey,  while  another 
[52]  discussed  the  identification  of  DHTs,  without  providing  an  implcmenation,  or 
suggesting  the  “triple”  heuristic  we  discuss  in  Section  3.2.2.  All  these  studies  required 
post-processing  on  the  graph  of  host  connections. 

Other  work  of  interest  includes  the  identification  of  compromised  “stepping 
stone”  computers  [56]  and  the  identification  of  relay  nodes  in  Skype  [51],  because 
such  hosts  are  qualitatively  similar  to  the  hosts  in  the  middle  of  the  “triples”  we 
define  in  Section  3.2.2.  The  identification  of  email  spamming  machines  by  graph 
theory  techniques  [11]  may  also  be  of  interest. 

One  might  think  of  a  “stepping  stone”  [56]  as  a  host  in  the  middle  of  a  triple 
(see  Section  3.2.2),  but  stepping  stone  analysis  requires  the  flows  to  be  very  similar 
on  the  packet  level.  In  the  interest  of  general  applicability,  we  only  concern  ourselves 
with  the  existence  of  the  flows,  no  with  any  of  their  measurable  properties  (other 
than  perhaps  their  overall  size). 


3.1.2  Trace  data 


We  use  publicly-available  traces  from  the  NLANR  Passive  Measurement  and  Anal¬ 
ysis  Project  [38].  The  PMA  project  provides  packet  header  traces  from  a  variety  of 
monitoring  locations.  Most  traces  have  a  duration  of  90  seconds,  but  special  traces 
with  durations  of  hours  or  days  are  also  available. 

Our  main  focus  was  on  the  Leipzig-I  data  set  [39],  a  trace  capturing  all  traffic 
of  the  University  of  Leipzig  between  8pm  Thursday,  November  21  and  2pm  Tuesday, 
November  26,  2002.  This  data  was  captured  recently  enough  that  peer-to-peer  traffic 
is  common,  bnt  sufficiently  long  ago  that  many  peers  were  still  using  fixed  TCP  and 
UDP  ports  [34],  In  addition  to  these  advantages,  the  long  time  period  of  the  trace 
allows  us  to  study  the  behavior  of  a  host  in  more  detail. 

We  can  also  compare  to  90-second  traces  from  the  University  of  Memphis  [40], 
captured  as  recently  as  March  2006. 

Because  they  each  monitor  the  gateway  to  a  university,  both  of  these  traces  have 
a  definite  “inside”  and  “outside.”  For  Leipzig-I,  the  inbound  and  outbound  packets 
are  stored  in  separate  hies  (whose  labels  were  inadvertently  switched  [33]).  We  place 
more  focus  on  classifying  the  interior  hosts,  since  the  trace  capture  more  information 
about  their  communications.  The  interior  hosts  may  also  be  of  more  interest  to  a 
network  administrator. 

We  may  look  at  what  happens  when  we  use  classifications  of  the  exterior  hosts 
to  adjust  our  judgements  of  the  interior  hosts. 
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3.1.3  Overview 


The  remainder  of  this  chapter  is  structured  as  follows.  In  Section  3.2  we  discuss 
the  measurable  features  of  our  data  set,  and  derive  a  model  for  message  passing.  In 
Section  3.3  we  discuss  the  use  of  port  numbers,  provide  an  overview  of  supervised 
learning  techniques,  and  list  the  results  of  our  classification  experiments. 
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3.2  Feature  set 


The  trace  data  we  use  provides  packet-level  records,  but  we  convert  them  to  NetFlow- 
style  data  [8]  befpre  performing  our  analyses.  In  particular,  we  use  the  following 
definition. 

Definition  3.1.  A  flow  is  a  sequence  of  packets  sent  back  and  forth  between  a  pair 
of  hosts,  identified  by  the  5-tnple 

/  —  (S,  s;  D,  d\  P),  (3.1) 

where  (S,  s)  are  the  source’s  IP  address  and  port  number,  (D,  d )  are  the  destination’s 
IP  address  and  port  number,  and  P  is  the  protocol  number.  The  “source”  is  the  host 
that  sends  the  first  observed  packet  in  the  flow.  If  there  is  a  gap  of  more  than  60 
seconds  between  consecutive  packets,  we  consider  there  to  be  two  distinct  flows. 

In  this  work,  we  only  consider  flows  corresponding  to  the  TCP  and  UDP  pro¬ 
tocols.  Most  other  protocols  could  easily  be  included,  but  for  the  triples  discussed  in 
Section  3.2.2,  it  is  important  to  ignore  ICMP  flows. 

Using  Definition  3.1,  we  divide  packet  traces  into  flows,  keeping  the  following 
information  about  each  flow: 

•  The  5-tuple  (S,  s;  D ,  d:  P),  with  s  =  d  =  0  for  protocols  without  port  numbers 

•  The  timestamp  of  the  first  packet  in  the  flow 

•  The  timestamp  of  the  last  packet  in  the  flow 

•  The  number  of  packets  sent  from  S  to  D 
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•  The  number  of  packets  sent  from  D  to  S 

•  The  number  of  bytes  (header  and  payload)  sent  from  S  to  D 

•  The  number  of  bytes  (header  and  payload)  sent  from  D  to  S 

•  A  flag  identifying  which  IP  address  (S  or  D)  is  internal  to  the  university 

3.2.1  Elementary  features 

Whenever  a  new  host  X  is  observed,  we  create  a  database  entry  for  it,  noting  the 
timestamp  of  the  first  observation  and  whether  the  host  is  interior  or  exterior  to 
the  university  (if  known).  As  X  participates  in  more  flows,  we  update  the  following 
values:  the  most  recent  timestamp;  a  list  of  the  distict  IP  addresses  X  interacts 
with;  the  total  number  of  flows  initiated  and  accepted;  the  total  number  of  packets 
initiated  and  accepted;  the  total  number  of  payload  bytes  initiated  and  accepted;  the 
total  duration  of  individual  flows;  and  counts  of  each  protocol  type. 

From  these  values,  we  calculate  nine  acceptable  features: 

•  The  average  number  of  new  neighbors  (IP  addresses  X  interacts  with)  per  sec¬ 
ond 

•  The  average  number  of  flows  per  neighbor 

•  The  average  number  of  packets  per  flow 

•  The  average  number  of  payload  bytes  per  packet 

•  The  average  flow  duration 
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•  The  fraction  of  flows  that  were  initiated  by  X 


•  The  fraction  of  packets  that  were  initiated  by  X 

•  The  fraction  of  payload  bytes  that  were  initiated  by  X 

•  The  fraction  of  flows  that  used  the  TCP  protocol 

These  features  can  be  calculated  at  a  specified  clock  time,  or  after  seeing  a  specified 
amount  of  activity  from  host  X. 

3.2.2  Triples 

We  now  present  a  feature  that  measures  flow  interactions. 

3. 2. 2.1  Definition 

For  a  decentralized  peer-to-peer  network  to  function,  the  peers  must  participate  in 
message-passing.  Whether  in  older  systems  such  as  Gnutella  or  newer  systems  based 
on  distributed  hash  tables  [32],  the  lack  of  a  central  authority  implies  that  when  a 
peer  wants  a  hie,  he  must  make  a  lookup  query  of  one  or  more  of  his  neighbors,  who 
must  in  turn  make  queries  of  one  or  more  of  their  neighbors,  until  the  hie  is  located. 

Therefore,  peers  in  a  decentralized  peer-to-peer  network  will  both  initiate  and 
accept  hows  with  other  peers,  and  in  quick  succession  if  they  are  not  willing  to  intro¬ 
duce  significant  delay  into  the  search  process.  This  inspires  the  following  definition, 
where  Time(/)  is  the  timestamp  of  the  hrst  packet  in  how  /. 


63 


Definition  3.2.  A  triple  for  a  host  X  is  a  pair  of  flows  fi  =  (A,  a;  X,  xp,pi),  f2  = 
(X ,  x2;  B ,  b;  p2)  such  that  A  /  B  and  0  <  Time(/2)  —  Time(/i)  <  r.  The  triple  is 
denoted  A  — >  X  — >  B.  The  parameter  r  is  called  the  window  parameter. 

We  use  a  window  parameter  of  r  =  5  seconds. 

Note  that  either  flow  f\  or  f2  may  be  involved  in  more  than  one  triple.  In  fact, 
in  the  case  of  Gnutella,  we  would  expect  f\  to  be  in  multiple  triples,  since  Gnutella 
lookup  queries  are  passed  to  all  neighbors. 

3. 2. 2. 2  Counting  triples 

Because  we  only  observe  traffic  that  takes  place  between  a  host  inside  the  university 
and  a  host  outside  the  university,  we  will  be  artificially  unlikely  to  see  triples  for 
exterior  hosts.  We  therefore  focus  on  interior  hosts.  We  will  find  that  P2P  hosts  tend 
to  have  more  triples,  even  though  we  cannot  see  the  peer  activity  that  takes  place 
entirely  within  the  university  network.  In  the  Leipzig-I  dataset,  we  can  easily  limit 
our  analysis  to  interior  hosts. 

Consider  observing  a  interior  host  X  for  a  time  period  0  <  t  <  T .  We  keep 
track  of  the  following  statistics: 

•  N\  —  the  number  of  flows  X  accepts 

•  N2  =  the  number  of  flows  X  initiates 

•  Si  —  the  number  of  accepted  flows  that  are  involved  in  a  triple 

•  S2  —  the  number  of  initiated  flows  that  are  involved  in  a  triple 
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H  =  the  number  of  distinct  hosts  X  communicates  with  (may  be  omitted 
see  below) 


The  flows  counted  by  N\  have  the  potential  to  be  the  first  flow  in  a  triple, 
while  those  counted  by  N2  have  the  potential  to  be  the  second.  We  will  model  each 
flow’s  probability  of  being  in  a  triple  by  assuming  that  the  start  times  of  the  flows 
are  uniformly  distributed  in  the  interval  [0,  T],  and  then  compare  this  expectation  to 
the  observations  S). 

We  choose  the  uniform  distribution  for  simplicity,  and  will  find  it  to  be  effective. 
We  tried  the  Poisson  distribution  with  similar  results,  but  it  has  been  known  for  years 
[43]  that  Poisson  models  do  not  capture  the  burstiness  in  data  networks.  Non-Poisson 
arrival  models  [41]  could  also  be  considered. 

Let  /1  be  one  of  the  N\  inbound  flows,  and  let  t\  =  Time(/i).  We  define 


Pi  =  P[/!  is  in  a  triple] 


=  P[a  flow  f2  to  a  different  host  occurs  with  Time(/2)  €  +  r] 


T\MW1n2 

T, 


(3.2) 

(3.3) 

(3.4) 


where  ^=lN2  is  the  average  number  of  flows  that  X  initiates  to  hosts,  other  than 
the  host  involved  in  f\ .  Similarly,  given  an  outbound  flow  f2,  we  can  define  the 
probability 

P2  =  1-(i-^)tA,‘  (3.5) 

that  f2  is  involved  in  a  triple. 
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If  H  is  deemed  to  be  to  expensive  to  keep  track  of,  we  can  use  the  approximations 

/  t  \N2  /  t  \  Nl 

pi«l-(l--j  ,  p2«l-(l--J  (3.6) 

because  ~  1  for  large  H. 

For  i  =  1,2,  we  can  now  model  the  number  of  flows  that  are  involved  in  a 
triple  using  a  Binomial(lVj,pj)  distribution,  if  we  assume  each  flow  has  an  indepen¬ 
dent,  identically  distributed  probability  of  being  involved  in  a  triple.  Recall  that  the 
binomial  distribution,  which  has  probability  mass  function 

Bi(n)  =  <3-7) 

has  mean  pt  =  NiPi  and  variance  of  =  Nlpl ( 1  —  pi). 

The  Binomial (A^j,  pj)  distribution  can  be  approximated  by  a  Normal (/i,,  of)  dis¬ 
tribution.  Therefore,  to  compare  5)  to  the  expected  /y ,  we  may  use  the  z-score 


S%  pi 


(3.8) 


The  z-score  represents  the  difference  between  S)  and  the  model  mean,  in  units  of  the 
model  standard  deviation.  The  ^-scores  z\  and  z 2  will  be  our  metrics  for  measuring 
triples. 
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3.3  Classification  results 


This  section  describes  our  testing  methodology  and  results. 


3.3.1  Determining  the  “true”  classes 

Although  TCP/UDP  port  numbers  are  not  generally  an  effective  method  for  identi¬ 
fying  peer-to-peer  traffic  [34],  we  can  use  them  as  a  baseline  measurement  of  truth 
in  the  Leipzig-I  data  set  because  of  its  age.  Table  3.1  summarizes  how  much  P2P 
traffic  can  be  identified  by  port  number;  we  see  that  port  numbers  cannot  be  used  as 
a  reliable  means  of  identification  in  the  Memphis  data  set. 


Table  3.1:  The  fraction  of  TCP/UDP  traffic,  by  number  of  ffows  and  by 
amount  of  payload,  that  has  known  P2P  port  numbers. 


Application 

Leipzig-Ia 
Flows  Bytes 

Memphis6 
Flows  Bytes 

eDonkey  2000 

23.54% 

12.38% 

0.28% 

0.01% 

Soribada 

1.52% 

0.04% 

— 

— 

WinMX 

1.19% 

3.58% 

— 

— 

DirectConnect 

1.05% 

1.15% 

— 

— 

KaZaA 

1.04% 

10.59% 

0.04% 

0.01% 

Gnutella 

0.48% 

1.27% 

0.75% 

0.07% 

Manolito 

0.21% 

0.01% 

0.73% 

0.01% 

SoulSeek 

0.10% 

0.65% 

— 

— 

Napster 

0.04% 

0.11% 

0.08% 

0.00% 

iMesh 

0.01% 

0.01% 

— 

— 

Bit  Torrent 

0.00% 

0.00% 

0.10% 

0.00% 

Client/Server 

Unknown 

48.51% 

22.29% 

35.21% 

34.99% 

50.94% 

47.08% 

54.34% 

45.55% 

“From  20021121-200000. 
bFrom  MEM-1143763417-1. 


Given  an  interior  host  X,  we  can  assign  one  of  four  labels  to  each  of  its  flows: 
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1.  Peer-to-Peer:  At  least  one  port  is  known  to  be  associated  with  a  peer-to-peer 
application  (KaZaA,  eDonkey  2000,  etc.) 

2.  Server:  At  least  one  port  associates  the  host  with  the  server  end  of  a  specific 
non-P2P  application  (HTTP,  games,  instant  messaging,  trojans,  etc.). 

3.  Client:  At  least  one  port  associates  the  host  with  the  client  end  of  a  specific 
non-P2P  application. 

4.  UNKNOWN:  Neither  port  has  a  known  association. 

We  use  the  port  list  in  Appendix  A.l,  derived  from  a  variety  of  online  sources. 
For  example,  if  X  uses  port  80  (HTTP),  we  label  the  flow  Server;  if  the  other  host 
uses  port  80,  we  label  the  flow  Client. 

If  at  least  75%  of  the  flows  that  host  X  participates  in  are  Unknown,  then 
we  also  assign  A"  the  label  Unknown.  Otherwise,  we  assign  X  one  of  the  labels 
Peer-to-Peer,  Server,  or  Client  by  a  simple  plurality  vote  of  the  relevant  flows. 
Counts  are  shown  in  Table  3.2.  Some  time-of-day  trends  are  apparent. 

We  conjecture  that  many  of  the  Unknown  flows  are  peer-to-peer,  because 
traditional  client/server  applications  tend  to  use  well-known  ports.  We  also  conjecture 
that  some  flows  taking  place  over  port  80  are  also  P2P  [26],  since  hosts  may  try  to 
hide  their  traffic  to  bypass  firewalls  or  other  security  measures. 


Table  3.2:  The  number  of  active  hosts  of  each  type,  broken  down  by  date  and 
time.  An  “active”  host  is  one  that  both  sent  and  received  one  packet  in  the  six- 
hour  file  window.  Not  all  hosts  were  used  in  the  experiments  (Section  3.3.3), 
because  the  experiments  employed  a  more  limited  time  (or  traffic  volume) 
window. 


Time 

Date 

P2P 

Client 

Server 

Unk. 

Total 

8pm 

Thu 

189 

2129 

312 

1120 

3750 

Fri 

144 

1492 

316 

513 

2465 

Sat 

133 

1205 

312 

301 

1951 

Sun 

190 

1880 

363 

344 

2777 

Mon 

226 

2364 

324 

454 

3368 

2am 

Fri 

102 

788 

169 

246 

1305 

Sat 

100 

415 

256 

309 

1080 

Sun 

104 

646 

236 

258 

1244 

Mon 

93 

768 

780 

98 

1739 

Tue 

122 

872 

203 

398 

1595 

8am 

Fri 

167 

3462 

388 

463 

4480 

Sat 

123 

1482 

214 

419 

2238 

Sun 

131 

1254 

857 

468 

2710 

Mon 

172 

3929 

561 

453 

5115 

Tue 

171 

3917 

520 

572 

5180 

2pm 

Fri 

184 

3103 

388 

432 

4107 

Sat 

139 

1683 

400 

352 

2574 

Sun 

169 

1866 

439 

220 

2694 

Mon 

227 

4326 

437 

535 

5525 

3.3.2  Supervised  learning  techniques 


We  will  compare  the  results  of  a  variety  of  supervised  learning  techniques.  This 
section  provides  an  overview  of  such  techniques.  Since  we  are  only  using  the  tools, 
not  developing  new  techniques,  we  do  not  provide  full  details;  see,  e.g.,  [24]  for  more 
information. 
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3. 3. 2.1  Goals  and  definitions 


The  goal  of  supervised  learning  is  to  construct  a  function,  based  on  known  examples, 
that  relates  observable  features  to  an  outcome  of  interest.  It  can  be  used  to  answer 
questions  such  as:  What  is  the  probability  of  a  customer  defaulting  on  their  mortgage? 
Will  a  patient’s  cancer  recur?  Is  a  computer  access  malicious? 

Each  event  (item  in  the  sample  space)  has: 

•  a  feature  vector  ~x  =  (aq,  x2 xp)  G  of  observable  variables1 

•  an  outcome  variable  y  G  Y2 

There  is  a  body  of  known  training  data  (aq,  yi), . . . ,  (xn,  y n).  Supervised  learn¬ 
ing  techniques  attempt  to  build  a  function  /  :  — >  Y  that  is  close  to  having 

f(xj)  =  yj  for  all  1  <  j  <  N,  without  overfitting.  Building  /  may  be  time-consuming, 
but  once  the  setup  process  is  complete,  predicting  the  outcome  y  =  f(x)  for  a  new 
observation  x  is  fast. 

3. 3. 2. 2  Examples 

We  begin  with  a  very  simple  example.  Assume  there  are  only  two  possible  outcomes, 

y  —  0  or  y  —  1.  The  training  data  divides  into  two  classes,  each  with  Ny  members, 
1It  is  possible  to  have  categorical  variables  £  Xt,  where  Xt  is  not  R,  nor  even  a  metric  space;  for 
example,  Xt  =  {female,  male}.  Since  none  of  our  features  are  categorical,  we  use  Xi  =  R  throughout 
to  simplify  the  discussion. 

2If  Y  is  finite,  building  /  is  known  as  a  “classification  problem”;  if  Y  =  R,  it  is  known  as  a 
“regression  problem.”  We  focus  on  classification  problems  here. 
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and  each  with  mean 


=  (3-9) 

V  3= 1 

and  covariance  matrix 

1  Ny 

Ey  =  w^ -ih)T(xj(y)  ~ih)-  (3-10) 

V  3=1 

If  the  true  classes  are  normally  distributed  with  equal  covariances,  then  Fisher’s  linear 
discrminant  [15]  provides  the  Bayes  optimal  solution.  This  discriminant  divides  the 
feature  space  M.p  with  a  hyperplane  perpendicular  to 


Tv  —  ((No  —  1)£0  +  (Ni  —  l)£i)  1(JTo  ~  hi)-  (3.11) 


That  is,  either 


fc(~x)  =  < 


if  w  ■  TE  >  c 
if  w  ■  x  <  c, 


or 


fc{~x ) 


I  1  if  w  ■  x  >  c 

0  if  Tv  ■  ~x  <  c, 


(3.12) 


for  some  constant  c. 

We  now  briefly  describe  some  of  the  most  important  supervised  learning  tech¬ 
niques.  In  a  given  problem,  there  may  be  more  than  two  possible  outcomes;  these 
techniques  can  be  generalized  to  multi-class  versions  with  varying  amounts  of  diffi¬ 
culty. 


•  Linear  discriminant  analysis  [15].  The  feature  space  is  divided  by  a  hyperplane. 
The  model  assumes  the  classes  have  equal  covariances. 


•  Quadratic  discriminant  analysis  [13].  The  feature  space  is  divided  by  a  quadric 
surface.  The  model  assumes  the  classes  are  normally  distributed  with  unequal 
covariances. 
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•  Naive  Bayes  classification  [12].  All  features  are  assumed  to  be  independent. 

•  k- nearest  neighbor  algorithm  [13].  A  new  observation  x  is  assigned  the  most 

common  y  value  of  its  k  nearest  neighbors  in  the  training  set. 

•  Neural  networks  [30].  The  function  /  is  a  compilaton  of  other  functions,  which 
fit  a  given  model  and  are  tuned  by  the  training  data. 

•  Classification  trees  [6].  The  feature  space  is  hierarchically  divided  along  one 
variable  at  a  time.  At  each  point  in  the  hierarchy,  the  split  is  chosen  to  reduce 
the  node  impurity  i  = 

•  Random  forests  [4],  A  large  number  of  classification  trees  are  produced,  each 
using  a  small  random  subset  of  the  features  and  a  randomly  chosen  (with  re¬ 
placement)  subset  of  the  training  data. 

We  use  Breiman’s  implementation  for  random  forests  [5],  and  LNKnet’s  imple¬ 
mentation  for  all  other  techniques  [23].  In  LNKnet,  we  performed  a  simple  mean-0, 
variance- 1  normalization  of  the  data;  such  normalization  is  unnecessary  in  random 
forests. 

3.3.3  Results 

Since  the  makeup  of  a  university’s  active  computer  population  varies  over  the  course 
of  a  day,  we  grouped  the  input  files  by  local  start  time:  8pm,  2am,  8am,  and  2pm. 
For  n  =  4, 16,  64,  256, 1024,  we  observed  the  first  n  flows  for  each  IP  address.  Again 
for  n  =  4, 16,  64,  256, 1024,  we  observed  the  first  n  seconds  of  activity  for  each  IP 
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address.  The  features  described  in  Section  3.2  were  computed,  then  run  through  var¬ 
ious  supervised  learning  algorithms.  To  account  for  imbalances  in  class  distribution, 
Peer-to-Peer  hosts  were  given  a  weight  five  times  as  great  as  Client  and  Server 
hosts  in  the  training  process. 

We  estimated  the  success  rate  (percent  of  Peer-to-Peer  hosts  classified  cor¬ 
rectly)  and  false  positive  rate  (percent  of  non- Peer- to- Peer  (Client  and  Server) 
hosts  classified  as  Peer-to-Peer)  for  each  algorithm.  With  random  forests,  these 
come  from  the  out  of  bag  (OOB)  error  estimate.  With  the  other  algorithms,  these 
come  from  10-fold  cross-validation  (CV10). 

A  comparison  of  classification  techniques  is  given  in  Tables  3.3  and  3.4.  We 
see  that  the  random  forests  technique  generally  has  the  highest  success  rates.  It  also 
gives  larger  false  positive  rates,  but  as  discussed  above,  we  believe  many  of  these 
“false  positives”  are  true  Peer-to-Peer  hosts  disguising  their  traffic  as  Client, 
Server,  or  Unknown  traffic. 

The  dependence  on  n  (the  number  of  flows  or  seconds)  is  depicted  in  Figure  3.1, 
for  random  forests.  We  see  that  increasing  the  number  of  flows  in  the  sample  increases 
the  success  rate  (and  the  false  positive  rate),  while  increasing  the  time  scale  of  the 
sample  has  little  effect. 

3.3.4  Testing  against  the  Memphis  data  set 

We  compare  Memphis  data  collected  just  after  midnight  on  Saturday,  March  31, 
2006  with  a  Leipzig  trace  collected  at  2am  on  Saturday,  November  23,  2002.  Because 
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Percent  Percent 


RandomForest  at  02:00 


(a)  Based  on  n  flows,  at  2am 
RandomForest  at  02:00 


(c)  Based  on  n  seconds  of  data,  at  2am 


RandomForest  at  1 4:00 


(b)  Based  on  n  flows,  at  2pm 
RandomForest  at  1 4:00 


(d)  Based  on  n  seconds  of  data,  at  2pm 


Figure  3.1:  Sample  of  the  results  for  the  random  forest  classifier.  Each  line 
represents  a  different  date.  The  upper  (black)  lines  are  success  rates.  The 
lower  (blue)  lines  are  false  positive  rates.  The  solid  lines  incorporate  all  fea¬ 
tures.  The  dashed  lines  incorporate  only  z\,  ,  and  the  fraction  of  traffic 

initiated  by  the  given  host. 
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Table  3.3:  The  range  of  success  and  false  positive  rates,  across  different  days, 
for  a  variety  of  classifiers,  based  on  the  first  16  flows  from  each  host.  The 
first  line  under  each  classifier  is  for  the  full  feature  set.  The  second  line  is  for 
the  limited  feature  set  consisting  of  z i,  Z2,  and  the  inbound/outbound  traffic 
ratios. 


Technique 

8pm 

S%  FP% 

2am 

S%  FP% 

8am 

S%  FP% 

2pm 

S%  FP% 

Random  forest 

62-75 

9-12 

77-81 

11-20 

48-65 

4-12 

56-65 

5-8 

50-67 

13-15 

71-79 

14-29 

41-55 

4-13 

43-57 

7-11 

Neural 

47-64 

6-12 

37-82 

2-40 

22-45 

1-9 

35-49 

3-7 

network 

37-54 

16-19 

44-78 

18-55 

15-37 

3-17 

16-33 

4-11 

Classification 

45-56 

5-7 

58-66 

5-12 

27-43 

2-7 

32-46 

4-5 

tree 

30-37 

7-8 

47-56 

8-18 

16-34 

3-8 

26-28 

4-8 

Linear 

42-52 

8-9 

37-69 

4-23 

8-34 

0-5 

24-43 

2-7 

discriminant 

33-42 

10-12 

30-59 

7-33 

1-27 

0-6 

13-28 

4-7 

Nearest 

35-40 

5-7 

45-64 

4-13 

22-34 

2-5 

26-41 

3-5 

neighbor 

19-31 

7-9 

30-53 

7-18 

13-20 

2-7 

15-21 

4-7 

Naive  Bayes 

27-41 

4-4 

24-57 

4-8 

16-29 

1-3 

23-35 

2-4 

2-12 

1-3 

5-39 

2-6 

0-13 

0-2 

2-11 

0-2 

port  number  are  not  a  reliable  means  of  identification  in  the  Memphis  data  set  (see 
Section  3.3.1),  we  use  Leipzig  as  the  training  data. 

Specifically,  we  build  a  random  forest  on  the  Leipzig  data  just  as  in  Section  3.3.3, 
then  pass  the  Memphis  data  through  the  forest.  We  base  the  random  forest  on  the 
first  64  seconds  of  data  from  each  host.  Less  than  2%  of  the  Memphis  hosts  had  16 
or  more  flows,  so  we  skip  that  comparative  analysis. 

Using  the  full  feature  set,  3%  of  the  Memphis  hosts  were  classified  as  Peer- 
to-Peer,  53%  as  Client,  and  43%  as  Server. 

Using  the  limited  feature  set  consisting  of  z±,  Z2,  and  the  inbound/outbound 
traffic  ratios,  which  we  consider  to  be  less  application-specific,  15%  of  the  Memphis 
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Table  3.4:  The  range  of  success  and  false  positive  rates,  across  different  days, 
for  a  variety  of  classifiers,  based  on  the  first  64  seconds  of  data  from  each  host. 
The  first  line  under  each  classifier  is  for  the  full  feature  set.  The  second  line 
is  for  the  limited  feature  set  consisting  of  z i,  22,  and  the  inbound/outbound 
traffic  ratios. 


Technique 

8pm 

S%  FP% 

2am 

S%  FP% 

8am 

S%  FP% 

2pm 

S%  FP% 

Random  forest 

67-78 

11-13 

81-89 

14-25 

64-79 

6-16 

69-72 

8-11 

44-65 

11-13 

66-84 

14-40 

51-75 

7-24 

49-62 

8-12 

Neural 

51-59 

3-10 

73-83 

12-56 

49-68 

2-27 

53-60 

3-8 

network 

55-64 

25-31 

62-93 

16-60 

57-70 

10-27 

48-63 

16-27 

Classification 

38-51 

6-9 

56-71 

7-16 

41-59 

2-9 

45-50 

5-7 

tree 

26-31 

6-9 

43-63 

7-36 

32-56 

3-18 

26-34 

5-10 

Linear 

35-50 

7-15 

50-83 

15-42 

36-63 

4-22 

34-43 

6-9 

discriminant 

34-47 

17-28 

47-88 

17-65 

35-49 

9-23 

20-48 

11-23 

Nearest 

29-53 

5-8 

51-68 

7-16 

38-54 

2-10 

30-45 

4-7 

neighbor 

23-30 

10-12 

39-59 

17-24 

27-37 

6-13 

25-33 

7-13 

Naive  Bayes 

18-39 

3-4 

37-53 

4-9 

36-56 

3-6 

33-38 

3-4 

12-19 

1-2 

28-49 

4-7 

22-36 

2-4 

14-21 

1-3 

hosts  were  classified  as  Peer-to-Peer,  71%  as  Client,  and  14%  as  Server. 

We  consider  these  results  promising.  The  two  data  sets  were  collected  more 
than  three  years  apart,  and  we  restrict  attention  to  only  interior  hosts  in  the  Leipzig 
data  set,  while  this  is  not  done  in  the  Memphis  data  set.  Nonetheless,  a  classifier 
built  from  the  Leipzig  hosts  is  sufficiently  sensitive  to  distinguish  different  behaviors 
in  the  Memphis  hosts,  rather  than  classifying  all  Memphis  hosts  as,  e.g.,  CLIENTS. 
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3.4  Discussion 


We  have  found  that  with  minimal  tweaking  of  the  parameters,  random  forests  can 
identify  80%  or  more  of  P2P  hosts.  One  study  [35]  added  a  variety  of  adjustments  to 
naive  Bayes  classifiers  to  improve  performance  from  65%  to  95%  on  certain  application 
categories;  if  they  had  started  with  a  more  powerful  classifier  such  as  random  forests, 
their  results  might  have  been  even  greater. 

Our  success  rates  are  comparable  to  those  of  related  studies,  given  the  limita¬ 
tions  on  our  data.  We  must  estimate  true  classes  from  port  numbers,  since  we  lack 
packet  payloads.  We  also  miss  all  peer  activity  taking  place  within  the  university;  we 
see  only  communcations  between  interior  and  exterior  hosts. 
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Chapter  4 


Discrete  analysis  of  crystal  steps 
4.1  Introduction 

In  Section  1.4,  we  discussed  the  importance  of  understanding  crystal  surface  evolution, 
particularly  the  motion  of  steps.  In  this  work,  we  develop  an  analytical  description 
of  step  motion  under  a  discrete  model  of  the  underlying  physics. 

Under  continuum  models,  the  crystal  surface  height  is  treated  as  a  continuous 
function  of  horizontal  position.  Because  physical  steps  have  discrete  heights,  discrete 
models  offer  a  more  natural  description  of  step  behavior.  Discrete  models  are  still 
approximations,  however,  because  the  horizontal  step  positions  are  typically  treated 
as  continuous  variables.  In  our  discrete  model  analysis,  we  find  an  extra  time  de¬ 
pendence  term  (f-1/2)  in  addition  to  the  exponential  dependence  that  emerges  in 
continuum  models. 

We  discuss  the  effect  of  the  model’s  physical  parameters  on  the  step  motion, 
and  therefore  the  step  bunching.  Consider  an  equilibrium  solution  corresponding 
to  equally  spaced  steps.  When  terrace  width  deviations  tend  to  decay  towards  this 
solution,  the  behavior  will  remain  stable.  When  the  terrace  widths  exhibit  growth 
away  from  equilibrium  instead,  step  bunching  can  occur. 
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4.1.1  Related  work 


A  discrete  model  for  crystal  steps  was  first  formulated  in  the  early  1950’s  [57],  after 
Burton,  Cabrera,  and  Frank  noted  the  growth  of  crystals  in  solutions  that  were  much 
less  supersaturated  than  demanded  by  contemporary  models.  They  concluded  that 
the  crystal  surface  holding  the  solution  must  not  be  smooth  at  the  nanoscale.  Such 
discrete  models  are  now  called  “BCF  models.” 

A  variety  of  extensions  to  the  original  BCF  model  have  been  offered.  Finite 
attachment-detachment  rates  at  step  edges  play  a  large  role  in  step  behavior  [66]. 
Step-step  interactions  based  on  the  surface  free  energy  were  considered  early  on  [63] ;  a 
variety  of  short-  and  long-range  interactions  have  been  considered  since  [67].  Modeling 
the  crystal  as  a  stack  of  concentric  circular  discs  allows  the  introduction  of  curvature 
and  step  line  tension  [76].  Analyses  of  more  general  shapes  have  been  limited  to 
single-layer  islands  [58]. 

An  electric  held  can  also  bias  the  terrace  diffusion  [80]  by  applying  a  direct 
current  to  the  adatoms,  which  have  a  net  charge  due  to  their  interactions  with  the 
crystal  surface.  The  inclusion  of  an  electric  held  in  a  BCF  model  [79,  77]  is  particularly 
important  because  application  of  a  direct  current  is  the  most  common  method  of 
inducing  step  bunching  experimentally. 

Step  bunching  [68]  can  appear  in  systems  limited  by  the  step  edge  attachment- 
detachment  rates  [66].  Experiments  show  that,  under  the  application  of  a  direct 
current,  behavior  of  the  bunching  instability  depends  in  a  complex  way  on  the  system’s 
temperature  and  on  the  direction  of  the  applied  current  [70,  82],  Step  bunching  can 
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also  be  induced  by  deposition  from  the  surrounding  vapor,  as  shown  experimentally 
[65]  and  numerically  [72],  Finally,  step  curvature  itself  can  also  induce  bunching  in 
the  model  [61,  62], 

In  this  work,  we  seek  to  understand  step  bunching  of  straight  steps  via  ana¬ 
lytical  examination  of  the  natural  discrete  equations,  rather  than  by  a  continuum 
approximation  [73]  or  numerically  [72,  71,  78,  79].  We  focus  on  the  importance  of  the 
applied  electric  field. 

Partial  differential  equations  have  been  used  to  study  step  behavior,  for  example 
scaling  laws  for  bunch  widths.  These  PDEs  can  be  derived  directly  from  continuum 
versions  of  surface  free  energy  and  chemical  potential  [75,  76],  or  they  can  be  derived 
as  the  continuum  limit  of  discrete  step  equations  [73].  The  PDEs  are  highly  nonlinear 
and  tend  to  be  of  fourth  order  (e.g.,  equation  (5)  of  [73]).  Most  analyses  have  employed 
linearized  versions  [74]  for  tractability.  The  linearization  is  performed  about  the 
solution  corresponding  to  equally  spaced  steps  (e.g.,  equation  (6)  of  [62]).  Growth  of 
step  bunches  due  to  electromigration  has  been  studied  in  a  nonlinear  context  (e.g., 
equation  (31)  of  [77]). 

4.1.2  Overview 

The  remainder  of  this  chapter  is  structured  as  follows.  In  Section  4.2,  we  construct  a 
discrete  model  for  the  physics  of  step  motion.  In  Section  4.3,  we  derive  coupled  equa¬ 
tions  for  terrace  widths  based  on  this  model.  In  Section  4.4,  we  find  the  linearization 
of  this  step  motion  equation.  Finally,  in  Section  4.5,  we  find  approximate  solutions 
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to  this  linearized  equation.  A  table  of  relevant  parameters  and  variables  appears  in 
Appendix  B.l. 
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4.2  Physical  models  of  step  motion 


In  this  section,  we  derive  a  discrete  model  for  step  motion  based  on  the  underlying 
physics.  We  will  analyze  the  behavior  of  a  set  of  coupled  linear  ODEs,  providing  a 
contrast  to  the  analyses  of  linear  PDEs  common  in  the  literature. 

4.2.1  Assumptions  for  the  discrete  model 

We  make  the  following  assumptions  about  the  geometry  of  our  model: 

(A3)  The  steps  are  straight  (Figure  1.2(b)) 

(A4)  The  step  sequence  is  monotone,  downwards  to  the  right 
(A5)  The  step  height  a  is  constant 

These  assumptions  permit  a  one-dimensional  model  of  the  crystal  surface,  with 
the  position  of  the  jth  step  edge  denoted  Xj(t).  The  jth  terrace  has  width  Xj+i(t)  — 
Xj(t),  or  normalized  width  Uj(t)  =  ^ ,  where  L  is  the  average  terrace  width. 

We  also  make  the  following  assumptions  about  the  kinetics. 

(A6)  Surface  diffusion  is  faster  than  attachment-detachment  (ADL  kinetics) 

(A7)  Step-step  interactions  have  energy  U  oc  1/ (terrace  width)2  (e.g.,  entropic 
or  elastic  interactions) 

(A8)  There  are  no  thermal  fluctuations  (low  temperature) 

Step  motion  proceeds  from  the  motion  of  adsorbed  atoms  (adatoms)  on  the 
crystal  surface.  See  Figure  4.1.  Let  Cj(x,t )  be  the  local  concentration  of  adatoms  on 
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Figure  4.1:  The  crystal  substrate  is  composed  of  bonded  atoms.  Adsorbed 
atoms  (adatoms)  interact  with  the  surface,  but  move  across  it.  If  an  adatom 
attaches  at  a  step  edge  ( xn ),  the  step  edge  advances  by  one  lattice  constant. 
Similarly,  an  adatom  can  detach  from  a  step  edge,  causing  the  edge  to  retreat. 


the  jth  terrace,  Xj  <  x  <  x]+\.  Also  let  Jj(x,t )  be  the  adatom  current  on  the  j 


th 


terrace. 


4.2.2  Discrete  terrace  equations 

The  model  we  derive  here  is  very  similar  to  the  straight-step  model  in  [62]. 

The  step  velocity  law  comes  from  mass  conservation  of  adatoms  near  the  step 


edge: 


dx3 

dt 


n 

a 


+  Jj~ l(xj  ) 


(4.1) 


where  is  the  atomic  volume  and  a  is  the  step  height. 


The  adatom  current  obeys  Fick’s  first  law  of  diffusion  and  is  also  affected  by 


electric  drift  with  speed  Vq: 


dc  • 

Jj(x)  =  -D-j£  +  v„Cj. 


(4.2) 


The  adatom  density  c,-  satisfies  Fick’s  second  law  and  is  also  affected  by  the 
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desorption  time  r  and  deposition  rate  R : 


dc^ 

dt 


d.L  1 

—  — Cj  R 


dx 

d2c- 
=  D—1 
dx 2 


r 


vo 


dx 


TC> 


R. 


We  will  assume  there  is  no  deposition:  R  =  0. 

For  tractability,  we  make  the  quasi-steady  approximation  [57] 


dci 

dt 


0. 


(4.3a) 

(4.3b) 


(4.4) 


This  approximation  is  not  rigorous,  but  it  agrees  with  physical  observations.  The 
justification  is  that  step  velocities  (which  depend  on  attachment-detachment)  are 
much  slower  than  terrace  adatom  diffusion,  so  the  concentration  Cj  relaxes  to  steady- 
state  between  step  adjustments. 


4.2.3  Discrete  step  boundary  conditions 

We  assume  linear  kinetics  at  both  the  left  and  right  edges  of  the  jth  terrace,  that  is, 

Jj ( xt )  =  ~k+{cj (tf )  -  CT)  ’  (4- ' 5a) 

Jj(x J+i)  =  k-  1)  “  cjh)  >  (4-5b) 

where  k+  and  k-  are  the  attachment-detachment  rate  coefficients.  If  k+,k-  are  un¬ 
equal,  we  say  there  is  an  “Ehrlich- Schwoebcl  barrier”  at  the  step  edge  [78]. 

The  equilibrium  adatom  density  is  given  by  the  Arrhenius  equation, 

C]q  =  c0e^/kBT.  (4.6) 

Here  /q-  is  the  step  chemical  potential,  which  is  the  change  in  surface  free  energy  when 
an  adatom  attaches  at  the  jth  step  edge.  Positive  values  of  /.q  indicate  a  propensity 
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for  a  step  edge  to  emit  adatoms,  while  negative  values  indicate  a  propensity  to  accept 
adatoms. 

Because  we  are  ignoring  step  curvature  (Assumption  A3),  the  chemical  potential 
depends  only  on  step-step  interactions.  Assuming  nearest-neighbor  effects,  we  have 

d 

H  (Xj,  xj+i)  +  U  (Xj_i,  Xj)\ ,  (4.7) 

where,  from  Assumption  A7, 

U(Xj,Xj+ 1)  OC  7 - - - ry.  (4.8) 

(Xj+1  -  Xj)2 

Then,  subsuming  the  constants  of  proportionality  into  g,  we  have 

j^9[(^y_(^_ 

kBT  \xj+ 1  -  xj )  \Xj  -  xj- 1 

r/  is  a  dimensionless  constant  that  describes  the  strength  of  the  nearest-neighbor  step- 
step  interactions. 
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4.3  Step  motion  equation  from  the  model 

Here  we  find  an  equation  for  ^  that  depends  explicitly  only  on  the  terrace  widths. 
We  will  see  that  when  steps  interact,  ^  depends  on  the  terraces  two  neighbors  away, 
but  when  step-step  interactions  are  eliminated,  ^  depends  only  on  the  terraces  one 
neighbor  away. 


4.3.1  Adatom  density  and  adatom  current 


Under  the  qnasi-steady  approximation,  the  PDE  in  equation  (4.3)  becomes  an  ODE: 


d 2c?-  dcj  1 

D  ~r\  2  ~V0  -  ~Cj  =  °. 

OXz  OX  T 


(4.10) 


which  has  solutions 


Cj{x)  =  B+eR+{x~xi)/L  +  BTeR-{x-x^/L 


(4.11) 


where  B ±  will  be  determined  by  boundary  conditions,  and 


R±  =—±  — ymy, 

2D  2D  v  ' 

V  =  - 

VqT 


(4.12) 

(4.13) 


are  dimensionless  constants. 


Substituting  equation  (4.11)  into  equation  (4.2)  gives 


JJx)  =  -9+B+efl+(r-"')/1  -  q-B-eR-{x-z’'l'L, 


(4,14) 


where 


D 


(4.15) 


We  will  find  it  convenient  to  set 


q^  =  q±  +  k-,  q±  =  q±-k+.  (4.16) 

Substituting  into  equation  (4.5)  from  equations  (4.14)  and  (4.11)  gives 


B+qi  +  B-qZ  =  -fc+Cf , 

B+q teR*“’  +  BJqteR-n‘  =  k-C]ll} 

where 

Xj-\- 1  Xj 
u>  =  — — 

is  the  normalized  width  of  step  j.  Now  we  can  solve  for  Bf,  hireling 

Bf  =  ±d(uj )  [k-q~C^.1  +  k+q+eR*UjC]q]  , 

where 

d H  =  +  -  jrr u  1  _  +  o  u- 
q+q_eH+u — q+qZeK~u 

Finally,  from  equation  (4.6),  we  have 

Cf  =  coe 9(u73-u7-0 

where  g  <C  1  is  the  strength  of  step  interactions. 


(4.17a) 

(4.17b) 


(4.18) 


(4.19) 


(4.20) 


(4.21) 


4.3.2  Terrace  width 


Via  equation  (4.1),  the  normalized  terrace  width  satisfies 


M,  = 


Xj+ 1  Xj 


n 

La 


(4.22a) 

(4.22b) 
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(4.22c) 


=  —  +q+Bf+1  +  q-Bj+1 

-  q+Bf(l  +  eR+Uj )  -  q-B~(l  + 

+  q+B+_ieR+u^  +  q_Bj_ieR-u^ 

=  +C]+2  [d(uj+1)k-(q+ql  -  g_g+)] 

+  C]l,  [d(uj+1)k+(q+q+eR-u^  -  q_qXeR+u^) 

-d(Uj)k4q+qZ(  1  +  e*^)  -  <z_g;(l  +  eR-u0)] 

+  Cf  [-d(Mj)A;+(g+gleR-^(l  +  eR+u0  -  g_g+eR+^(  1  +  eR-u0) 
+d(Mj_i)fc_(g+gIei?+^-1  - 
+  C^!  [d(^_i)fc+(g+^eR-^-1eR+^-1  -  g.g+e^'-^^-^-1))]  (4.22d) 

=  ~LC°  e^Uj+2~Uj^  [-d(uj+1)k+k.{q+  -  g_)] 

e9{ui+1~ui  )  [+d(uj+i)k+(—q-qXeR+Uj+1  +  q+qXeR~Uj+1) 

-d(Uj)k-(q+qZ(  1  +  e*^)  -  <z_g;(l  +  e*-“'))] 

+  Au73~u7-i)  [-d(Uj)k+(k-(q+  -  q-)e{R++R ~)u> 

-  q_qleR+Uj  +  q+qXeR~Uj) 

+d(uj_1)k_(q+qZeR+ui-1  -  q^q-eR-u^1)] 

+  e9(uj-i-uj-z)  [+d(uj -i)k+k_(q+  —  g_)e(i?++i?-)uJ-i]  .  (4.22e) 

Therefore,  ^  depends  on  the  two  terraces  to  either  side:  Uj ,  Wj+i,  and 


Uj+2,Uj-2-  However,  the  dependence  on  Wj+2  and  Uj- 2  is  only  through  the  step-step 
interactions.  If  g  —  0,  ^  depends  only  on  Uj+1,  Uj,  Uj-\. 


4.4  Linearization 


To  make  progress  with  the  solution  of  equation  4.22e,  we  find  the  linearization.  The 
coefficients  of  u]+2  and  Uj_ 2  will  be  O  (g),  reflecting  that  the  two-away  neighbors 
only  affect  the  step  motion  through  step-step  interactions.  We  then  investigate  the 
dependence  of  all  the  coefficients  on  the  physical  parameters,  in  the  limiting  cases  of 
very  strong  and  very  weak  electric  held. 


4.4.1  Subordinate  functions 

We  will  linearize  around  the  solution  corresponding  to  equally  spaced  steps,  that  is, 
Uj  =  1  for  all  j.  We  define  the  deviation 

Uj  =  Uj  —  1  (4.23) 


so  we  may  consider  Uj  <C  1.  Several  parts  of  equation  (4.22e)  must  be  linearized. 


u 


eKuj  ~  eK  +  KeKUj  +  O  (tf)  , 

—  'Uj—i  ~  —  3 Uj  +  (D  ( Uj  +  Hj_i)  , 

d(uj)  r^j  do  d\Uj  T  O  i 


(4.24) 

(4.25) 

(4.26) 


where 


dn  — 


q+q_eR+  — q+q_eR- 


,  d\  = 


— q\q-R+eR+  +  q+qhR^eR~ 
ZeR+  —  q+qte^Y 


We  note  that 


(4.27) 


doR+  +  di  =  —  do(R+  —  R-)q+q^eR~ , 


(4.28a) 


doR-  “I-  d\  — 


(4.28b) 


-d20(R+  ~  R-)qXq_eR+ , 
d0(R+  +  R- )  +  di  =  +d20(qXqZR-eR+  -  q+q±R+eR~) 


(4.28c) 


4.4.2  Terrace  width  evolution 

We  now  find  the  linearization  of  Uj  =  iij.  We  require  each  of  \u3\  <C  1,  | guj\  <C  1,  and 
\R±Uj\  -C  1  for  all  j. 


U-: 


n 

La 


c0  |  (1  -  3gUj+2  +  3guj+i)  [-(d0  +  d1uj+1)(k+k-(q+  -  g_))] 

+  (1  -  3gUj+\  +  3 guj)  [+(d0  +  d1uj+1)(k+(-q_qXeR+  +  q+qteR-) 
+  k+(-q~qXR+eR+  +  q+qtR-eR-)uj+1 ) 

-  (d0  +  d1uj)(k-(q+ql(  1  +  eR+)  -  <?_g+(l  +  eR~ )) 

+k-(q+qZR+eR+  -  q~q+R~eR-)uj)] 

+  (1  -  3 guj  +  3guj_l)  [-(d0  +  d1uj)(k+(k_(q+  -  q~)eR++R- 

—  q_qXeR+  +  q+q+eR~) 

+  k+(k_(R+  +  R-)(q+  -  q-)eR++R- 

-  q-qXR+eR+  +  q+q±R-eR~)uj ) 

+  (d0  +  diUj-i)(k-(q+qZeR+  -  q-q~eR~) 

+k-(q+qZR+eR+  -  q-q+R-eR~)uj-i)\ 

+  (1  -  3guj_i  +  3 gUj-2)  [+(do  +  dluj_1)(k+k_(q+  -  q_)eR++R~ 


n 

La 


c0 


+k+k_(R+  +  i?„)(g+  -  q-)eR++R~Uj_i 
(1  -  3gUj+2  +  3 guj+i) 


+  0(u2)  (4.29a) 
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[~d0k+k-(q+  -  q-) 

-  dik+k-(q+  -  q-)uj+i ] 

+  (1  -  3guj+i  +  2>guj) 

[+d0  (k+k_(q+  -  q_)  (l  +  eR+  +  eR-) 

~(k+  +  /c_)g+g_  (eR+  -  eR “)) 

+  d20k2+(R+  -  R-)(q+  -  q-)qXqteR++R~uj+l 
+  k-(q+  -  q-)  [dlk-(R+  -  R~)qZrqZeR++R-  +  d±k+ )  Hj] 
+  (1  -  3 guj  +  Zguj-i ) 

[~d0  (k+k_(q+  -  q_)  ( eR++R-  +  +  eR-) 

~(k+  +  /c_)g+g_  (eR+  -  eR~ )) 

-  M<?+  -  <?-)  {d20k+(R+  -  R~)qXq± 

+k-(d0(R+  +  R-)  +  di))  ei?++'R_-u:,- 

-  d20k2_(R+  -  i?-)(g+  -  q-)q+qZeR++R-Uj- J 


+  (1  -  3guj-i  +  3gUj-2) 

[+d0k+k_(q+  -  g_)eR++R- 


+  (d0(-R+  +  -R_)  +  di)A;+A;_(g+  -  q_)eR++R~Uj_ i 


+  0  («2) 


(4.29b) 


O 


uj+ 2  [+3gd0k+k-(q+ 


+  uj+1  [(<?+  -  q-)d20{k\(R+  -  R_)qXqteR++R~ 

+  k+k_(qXqZR+eR+  —  qZrqJZR~eR ")) 
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-  3 gd0(k+k-(q+  -  q~)(eR+  +  eR~  +  2) 

-  (k+ +  k-)q+q.(eR+ -  eR-))] 

+  Uj[(q+  -  q-)d20(-(R+  -  R^){k2+qXqt  ~  k2_q~qZ)eR++R- 

-  k+k-((qXqlR~eR+  —  q+qJZR+eR~)eR++R~ 

+  (qXqlR+eR+  -  qZrqtR-eR~))) 

+  3gd0(k+k-(q+  -  q~)(eR++R-  +  2eR+  +  2eR~  +  1) 

-  2 (k+  +  k-)q+q_(eR+  -  eR~ ))] 

+  u3-\  [(q+  -  q~)d20eR++R-  {-k2_(R+  -  R-)q+qZ 

+  k+k_(qXqlR~eR+  —  q+qlR+eR~)) 

—  3gd0(k+k-(q+  —  q_)(  2eR++R~  +  eR+  +  eR~) 

-  (. k+  +  k_)q+q_(eR+  -  eH~))] 

+  Uj-2[+3gd0k+k-(q+  -  q-)eR++R~]  +  O  ( u 2)  .  (4.29c) 

Note  that  all  affine  terms  have  cancelled  out. 

We  will  typically  use  the  more  compact  notation 

Uj  ~  (gfo)  Uj+ 2  +  («i  +  gfii)  uj+i  +  («o  +  g/3 o)  uj 

+  (oi~  i  +  <7/3- 1)  hj-i  +  (g(3-  2)  (4.30) 

where 

«i  =  ^(9+  -  </_) do (/;+(/?+  -  R„)qlqteR++R- 

+  k+k_(qXqZR+eR+  —  q+qlR~eR~)),  (4.31a) 

a-i  =  7^(?+  -  q-)d20eR++R~ (~kX(R+  -  R-)q+qZ 
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+  k+k-(qXq_R~eR+  -  q+q±R+eR~)), 


(4.31b) 


«o  =  —  («i  +  ck— i ) ;  (4.31c) 

02  =  3 ^-d0k+k-(q+  -  q_),  (4.31d) 

P-2  =  3^-dok+k-  (g+:  -  q-)eR++R- ,  (4.31e) 

La 

Or 

ft  =  -3-^<f0(ft-M9+  -  9-)(eR+  +  eR-  +  2) 

-(ft  +  ft)<,+g_(eR*-eR-)),  (4.31f) 

O/- 

/3-i  =  -3  —^-d0(k+k-(q+  -  <?_)(2eR++R-  +  eR+  +  e^”) 

La 

-  [k+  +  k-)q+q-{eR+ -  eR-)),  (4.31g) 

Po  =  —  (P2  +  Pi  +  P- 1  +  P-2)  ■  (4.31h) 


4.4.3  Dependence  on  the  physical  parameters 

Here  we  examine  the  behavior  of  each  of  the  coefficients  an,  Pn.  A  summary  of  the 
physical  parameters  appears  in  Appendix  B.l.  We  should  always  have 

(A9)  L  >  0,  a  >  0,  f!  >  0,  Co  >  0,  r  >  0,  D  >  0,  k+  >  0,  k-  >  0 

We  will  also  assume 

(A10)  uo  >  0,  g  >  0 

Then  we  have  q+  >  0  >  g_  and  R+  >  0  >  R-. 

From  these  assumptions,  we  first  note  that 

do  =  —  qpqPeR~]  1  (4.32a) 

=  [( q+q _  -  k+q+  +  /c_g_  -  k+k-)eR+ 
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(4.32b) 


-  ( q+q _  +  /c_g+  -  /c+g_  -  k+kJ)eR  ] 

=  [{q+q-  -  k+k-)  ( eR+  -  eR-) 

—  q+  (k+eR+  +  k_eR~)  +  g_  ( k_eR+  +  k+eR~)]  1  (4.32c) 

<  0.  (4.32d) 


Holding  all  the  parameters  L,  a,  12,  c0,  r,  D,  k+,  k_  constant,  and  letting  v0  vary 
in  (0,  cx)),  we  have 


(  ( 


Lv0/D 

for  v0  1 

i?_  ~  < 

- L/(v0t ) 

1 

for  v0  1 

(4.33a) 

1 

for  ti0  «  1, 

1  A 

1 

✓ 

for  t)0  «  1, 

<?+  ~  < 

( 

D/ (vqt  ) 

for  Vq  S>  1 

q-  ~  < 

/ 

-u0 

for  t’o  »  1 

(4.33b) 

/DA 

for  !)0  «  1, 

-VdF 

V 

for  j)0  «  1. 

Further  assuming  that  5  —  L/ \J Dt  -C  1, 

we  have 

d0  ~  S 

/ 

—  [k_v0eLv°/D]  1 

i 

for 

v0  >  1 

(4.34) 

2(k+  +  k-)y/D/r 


+  0(5) 


for  u0  «  1. 


\ 


4. 4. 3.1  Non-interaction  terms  (0(1)) 


From  equation  (4.31a),  we  can  calculate 


Oi\ 


Oc0 

La 


fc+(fc++fc_)  Lg 0_e-Lvo/D  for  ^  1 


■  < 


-8JH: 


fc_  D 

D  k+ 

T  k++k- 


(4.35) 


for  v0  <  1. 
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From  equation  (4.31b),  we  can  calculate 


( 


CM.— I 


VLc0 

La 


•  < 


k+~ 

'  VQ  T 


for  t’o  »  1 


(4.36) 


(6\/vkLLkI  iorv  o<l- 

When  the  electric  field  is  very  strong,  a±  — *  0”  and  a_i  — >  0+,  with  an  decaying 
faster  in  v0  than  a_i. 

In  the  case  where  there  is  no  Ehrlich-Schwoebel  effect,  that  is,  k+  =  k-,  then, 
when  the  electric  field  is  very  weak,  we  have  an  ~  — a_i,  and  a0  ~  0. 


4. 4. 3. 2  Step-step  interaction  terms  (O  (g)) 


From  equation  (4.31d),  we  know  /?2  <  0  for  all  valid  choices  of  the  physical  parameters. 
We  can  calculate 


a  Qc0  . 


—3 k+e  Lv°/D  for  t’0  1 


(4.37) 


-3 


k+k— 

k+-\-k— 


for  u0  <C  1. 


From  equation  (4.31e),  we  know  it  is  also  true  that  /3_2  <  0  for  all  valid  choices 
of  the  physical  parameters.  We  can  calculate 


0-2 


flco 

La 


-3 k+  for  t’o  S>  1 


-3wfc  for  u0  «  1. 


(4.38) 


From  equation  (4.31f),  we  know  [3\  >  0  for  all  valid  choices  of  the  physical 
parameters.  We  can  calculate 


n  tlCo  . 

01  Ta  { 


3k , 


for  t’o  3>  1 


(4.39) 


12WTW  for  u0  «  !. 
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From  equation  (4.31g),  we  know  /3_i  >  0  for  all  valid  choices  of  the  physical 


parameters  as  well.  We  can  calculate 


P- 1  ~ 


h2co 

La 


for  v0  S>  1 


(4.40) 


[12WTW  for  Uo  -C  1. 

When  the  electric  held  is  very  strong,  we  have  p2  ~  P-2,  while  /5_i  m  3/3i . 
When  the  electric  held  is  very  weak,  we  have  p2  ~  P-2,  and  Pi  ~  P-\. 
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4.5  Solutions  to  linearized  equation 

We  transform  the  linearized  step  motion  equation  to  simplify  the  coupling,  thus  find¬ 
ing  an  integral  expression  for  Uj(t).  We  then  use  the  method  of  steepest  descents  to 
estimate  this  integral  for  long  times  t. 

4.5.1  Transformation 

Emulating  the  solution  to  a  difference  equation  in  [59],  we  assume  that  Uj(t)  takes 
the  form 

Uj(t)  =  f  e^z  f  (z,  t)dz  (4.41) 

JC(t) 

for  all  j,  t,  for  some  fixed  function  f(z,  t )  and  contours  C(t).  Writing  the  endpoints  of 
C(t)  as  Co(t),  C'i(t),  and  assuming  we  may  interchange  integration  and  differentiation, 
the  Leibniz  integral  rule  gives 

fi j(i)=  [  (»),»)  -  £*(»)«<#*»/(£*(»),«).  (4.42) 

JC(t)  at 

To  eliminate  the  last  two  terms,  we  require  the  endpoint  contribution  of  the  integrand 
to  vanish,  so  e^zf{z,t)  — >  0  in  measure  as  z  — *  a(t),  z  — ■>  b(t). 

For  convenience,  we  define 

h(z)  =  {g^e21"  +  (cki  +  g(3i)elz  +  (cto  +  <?A))  +  («-i  +  9@-i)e  tzd-(gP- 2)e  2t~-  (4.43) 
Now,  substituting  equation  (4.41)  into  equation  (4.30)  gives 

[  eijz?f{z,t)dz~  [  eijzh{z)f(z,t)dz.  (4.44) 

JC(t)  Ot  Jc (t) 

For  this  equation  to  hold  for  all  j,  f(z,t )  should  satisfy  the  differential  equation 

^(z,t)  =  h(z)f(z,t),  (4.45) 
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which  has  solutions 


f(z,t)  =  U{z)eth{z\  (4.46) 

We  will  use  the  requirement  U{x)  — >  0  as  x  — >  ±oo  to  ensure  the  vanishing 
endpoint  contributions  in  equation  4.42.  Then  we  can  take  C(t)  =  M  for  all  t,  so 
that  equation  4.41  is  a  simple  Fourier  transform.  The  other  factors  el^z  and  eth ^  are 
oscillatory  along  this  contour,  confirming  decay  at  the  endpoints. 

4.5.2  Initial  data 

U(z)  can  be  determined  by  the  initial  data  of  the  problem.  We  have  Uj(0)  = 
fm  e^zf(z>  0)dz  =  fR  elizU(z)dz.  Treating  j  as  a  real  (rather  than  integer)  parameter, 
we  see  that  U(z)  is  the  inverse  Fourier  transform  of  Uj( 0): 

U(z)  =  4  f  e-“%(0 )dj.  (4.47) 

Jr 

To  represent  a  valid  physical  system,  Uj{ 0)  must  be  continuous  in  j,  the  average 
of  the  uj( 0)  over  all  j  &  Z  must  equal  0,  and  we  must  have  Uj( 0)  >  —  1  for  all  j.  We 
also  restrict  ourselves  to  {q(0)  such  that  U(z)  — >  0  as  z  — >  ±oo  in  a  neighborhood  of 
the  real  axis. 

Later,  we  will  need  to  use  the  fact  that  that  U(—z)  =  U{z),  which  is  easily 
verified: 


[  e^uMdj 

Jr 

(4.48a) 

1 

”  2tt  , 

1 

“  2 n. 

[  er^uj(0)dj 

Jr 

r 

(4.48b) 

/  e~vJzUj{d)dj 

Jr 

(4.48c) 

(4.48d) 


=  U{z) 

As  a  consequence,  U  maps  the  imaginary  axis  to  the  real  axis,  since  U (■ iy )  =  U (— iy )  = 
Wry). 

As  an  example,  we  might  choose  Uj( 0;  e)  =  which  yields  U(z]e)  = 

—  ^  •  Since  Uj(0]  e)  is  an  odd  function  of  j,  we  are  guaranteed  that 

1  N 

— I  £  «i(°;c)  =  0  (4.49) 

j=-N 

for  all  N,  including  N  — >  oo. 

4.5.3  Steepest  descents 

From  equations  (4.41)  and  (4.46),  we  have 

Uj(t)  =  [  U{z)eijzeth{z)dz,  (4.50) 

Jr 

where  U{z)  is  given  by  equation  (4.47).  In  each  case  below,  we  will  consider  the 
complex  plane  restricted  to  x  G  (— 7r,  7r]  to  determine  the  relevant  critical  point(s)  and 
contour  of  integration.  This  region  should  be  repeated  at  intervals  of  2mr  (n  e  Z)  to 
deform  the  entire  real  line  M.  e*JZ  and  eth<'z'1  are  periodic  in  Real  (z)  with  period  2ir, 
so  only  the  value  of  I4(z)  will  differ  between  regions. 

We  are  interested  in  the  long-term  behavior  of  the  terrace  widths,  so  we  assume 
t  1.  We  consider  steps  “in  the  bulk”  of  the  crystal,  where  |^|  <  1.  To  approximate 
the  above  integral  using  Riemann’s  method  of  steepest  descents,  we  must  locate  the 
critical  points  where  h\z)  =  0.  From  h(z)  in  equation  (4.43),  we  have 

h\z)  =  i  (2 (9p2)e2iz  +  (ai  +  g(3 i)eiz  -  (a_x  +  g(3-i)e~iz  -  2(gf3_2)e~2iz)  .  (4.51) 
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Letting  w  =  etz,  we  know  that  w  cannot  equal  0.  Therefore,  h'(z)  =  0  if  and  only  if 
w  is  a  nonzero  root  of  the  polynomial 

pg{w)  =  (2 gP2)w4  +  («i  +  g(5i)wz  -  (cm i  +  g/3-i)w  -  (2 (4.52) 


Suppose  we  have  determined  a  contour  C  passing  through  a  unique  maximum 
critical  point  z0,  where  h'(z0 )  =  0.  We  can  approximate  equation  (4.50)  by  making 
the  change  of  variable 


C2  =  h{zo)  ~  h(z) - ^h"(z0)(z  -  z0 )2  +  0((z-  z0 )3)  (4.53) 

so  that  z  ~  z0  +  \—\h"{zQ)\  1/2  C  +  O  (C2)  and  ^  ~  [-\h"{z0)\  1/2  +  C>  (C).  (Here 
we  must  £x  a  branch  of  the  square  root  function.  We  choose  it  such  that  \/^  =  \fz 
for  all  z  not  on  the  negative  real  axis.)  Then 


u 


j(t)  =  /  U(z)eijzeth{z)dz 


<c 


dz 


W(^(C))eiiz(?)et(/l("o)-<:2)^dC 


U(z0)eijz°eth{zo) 

1 


=  U(z0)  •  Vvr 


-/'(*)) 

-1/2 


-,-1/2 


e 


•  (eiz°)J  •  e^bo)/-1/2 


(4.54a) 

(4.54b) 

(4.54c) 

(4.54d) 


The  factor  is  a  departure  from  continuum  estimations.  It  is  the  fastest  decay 
we  could  expect  to  see  in  an  asymptotic  expansion.  If  the  first  non-zero  derivative  of 
h  were  of  higher  order  (h'(z0)  =  h"(z0)  =  •  •  •  =  h(m_1^(^0)  =  0,  h^m\z0)  ^  0,  m  >  2), 
the  decay  would  be  f_1//m. 
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(a)  Both  real 


(b)  Both  complex 


Figure  4.2:  The  case  g  =  0.  Contours  of  constant  imaginary  part  equal  to  the 
imaginary  part  of  a  critical  point,  shown  for  x  £  ( — 7r,  7t]  . 

4. 5. 3.1  No  step  interactions  ( g  =  0) 

In  this  case,  the  polynomial  (4.52)  is  simply 

Po(tv)  =  a\w 3  —  a-iw.  (4.55) 

We  must  ignore  the  root  w  =  0,  because  0  =  w  =  elz  has  no  solutions  for  finite  z. 
The  other  roots  are  ±Wq,  where  Wq  =  Depending  on  the  physical  parameters, 

Wq  may  be  real,  pure  imaginary,  or  zero.  Wq  =  0  is  a  bifurcation  point,  so  we  analyze 
the  other  two  cases  in  depth.  They  are  generically  pictured  in  Figure  4.2. 

The  deformed  integration  path  should  be  chosen  so  that  Imag  {h(z))  is  constant 
and  Real  (h(z))  reaches  a  maximum  at  a  critical  point.  We  have 

Imag  {h{x  +  iy))  =  sin(x)  ( a\e~y  —  a_iey)  ,  (4.56a) 

Real  (h(x  +  iy))  =  cos(x)  (aie~y  +  a_iey)  +  a0.  (4.56b) 


4. 5. 3. 1.1  Real  roots  (sign(ai)  =  sign(a_i)).  When  ta0  is  real,  h(zo)  is  real  as 
well,  so  we  must  find  the  contours  where  Imag(/r(^))  =  0.  We  immediately  see  that 


101 


Imag  (h(x  +  iy))  —  0  if  and  only  if  x  =  mi  (for  some  n  G  Z)  or  e~y  =  \J 

When  x  =  2n7r,  we  have  cos(x)  =  1,  so  Real  ( h(x  +  iy))  — >  —  oo  as  y  — >  ±cx). 
Similarly,  when  x  =  (2n  +  l)7r,  Real  (h(x  +  iy))  — >  +oo.  Therefore,  along  the  hori¬ 
zontal  contour  y  —  —  In  hr  Figure  4.2(a),  we  know  that  Real(/i(z))  reaches  a 

maximum  at  the  points  zn  =  (2 n  +  l)7r  —  i  In  We  have  z0  =  n  —  i  In 

For  the  final  calculation,  note  that 


eiza  _  _Wq  —  —  —L  <  o, 


h(z0)  =  +  Vl^-il)  >  0, 

h"(zo)  =  —  2^/aiQ!_i  <  0. 


(4.57a) 


(4.57b) 

(4.57c) 


Also  recall  that  U[—z)  =  U[z).  Now,  by  summing  equation  (4.54)  over  all  adjacent 
regions  of  width  27r,  we  have 


T  1  i 

(t)  ~  ^2u{z0  +  2wr)  •  Jn  -- h'\z0  +  2nn)  •  (ei{zo+2n7r))j  •  eM*o+2nn)t-i/2 


(4.58a) 


=  \/vr  —  ^/r"(^0)  •  (— wo)J  •  eth^z°H  l^2  ■  ^  7/  ^(2n  +  l)7r  —  i  In  \J~^j 


(4.58b) 


=  v71"  (aitt- 


^-i/4  _  ^  /Q-iy  .  et(v^T+Vl«-il)  ^ — 1/2 


2  Real  f  7/  f  (2 n  +  l)7r  —  i  In 


(4.58c) 


4. 5. 3. 1.2  Pure  imaginary  roots  (sign(ai)  ^  sign(o!_1)).  Here  w 0  =  i*  , 


so  =  (n  +  |)7r  —  i  In  *  1  for  n  e  Z  are  all  the  critical  points. 
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When  wq  is  imaginary,  h(zo)  is  complex.  Furthermore,  we  have  the  imagi¬ 
nary  part  Imag  (h(x  +  iy))  =  0  only  when  x  =  nn  (n  e  Z),  and  not  along  the 
horizontal  contour  in  Figure  4.2(a).  Because  the  contours  where  Imag  (h(x  +  iy))  = 
Imag  (h(z0)),  a  nonzero  constant,  may  not  cross  the  contours  where  Imag(/i(z))  =  0, 
they  are  restricted  to  regions  where  x  G  (m r,  (n  +  l)7r),  so  we  have  y  — >  ±oo.  There¬ 
fore,  to  keep  Imag  (h{x  +  iy))  from  growing,  we  must  have  sin(x)  — >  0  with  exponen¬ 
tial  decay.  This  implies  that  the  contours  are  asymptotic  to  x  —  nn  (n  e  Z). 

When  x  ~  2nir,  we  have  cos(x)  >  0,  so  Real  (h(x  +  iy))  — ■>  —  oo  as  y  — > 
— cx)  and  Real  (h(x  +  iy))  — >  +oo  as  y  — ■>  +oo.  Similarly,  when  a;  ~  (2n  +  l)7r, 
Real  (/i(x  +  iy))  — >  +cx)  as  y  — >  —  cx)  and  Real  (/i(x  +  iy))  — >  — cx  as  y  — >  +oo.  There¬ 
fore,  along  the  contour  joining  (— 7T,  +cx)  to  (0,  — cx)  to  (7r, +oo)  in  Figure  4.2(b), 
Real  ( h(z ))  reaches  maxima  at  the  critical  points  Zo,  —  2<j. 

For  the  final  calculation,  note  that 


eizo 


h(z0)  =  |oci 


(X—l 


2iV\ 


h"(z0)  =  2iy/|a1a_i 


(4.59a) 

(4.59b) 

(4.59c) 


Also  note  that  h(—z)  =  h(z)  as  well.  Now  we  have 


u 


nEZ 


l 


v)~£  U  (zq  +  2nn)  ■ 

^  W  (— zjj  —  2wr)  • 


+  2wr) 


-1/2 


nE  Z 


—  ^/i"(-z0-  -  2wr) 


(e*(z0+2n7r))/  _  gl/i(^0+2ri7r)^— 1/2 
1  -1/2 


(e*(-z0-2rwr)y  _  eth(-z0-2nn)^-l/2 


(4.60a) 
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7lEZ 


y/nt  1 y^(u(zp  +  2nn)  • 


—2h\z0) 


1-1/2 


,3  pth(z0) 


wOe 


U(zq  +  2nn)  ■ 


~2h"(zo) 


-1/2 


2y/nt  1//2  Real  f  W(z0  +  2n7r) 


nGZ 


waieth(zo) 


■ 2h"(z0 ) 


-1/2 


=  2v/vrt 


-1/2 


Real  (  i-\/|aiaZiT j  ': 


wjneth{zo) 


Jh(zo) 


■  Real  (U  (z0  +  2n7r)) 

nEZ 

I  mag  ^-^laia-il) 

•  Imag  {U  (z0  +  2mr)) 

nEZ 


a_i 


( y,\ 


r,th(zo) 


(4.60b) 

(4.60c) 


(4.60d) 


4. 5. 3. 2  Step  repulsions  ( g  >  0) 

In  this  case,  pg(w )  is  given  by  the  full  equation  (4.52).  Note  that  pg(w )  has  real 
coefficients,  so  if  Wq  is  a  root,  then  Wq  is  also  a  root.  This  implies  that  if  z0  is  a 
critical  point,  then  —Zo  is  also  a  critical  point. 

We  know  that  (3± 2  are  both  negative.  Therefore,  pg(w )  — >  —  00  as  w  — >  ±00, 
and  ps(0)  >  0.  The  intermediate  value  theorem  implies  that  pg(w)  has  at  least  two 
real  roots,  one  positive  and  one  negative.  Furthermore,  if  pg(w )  has  four  real  roots, 
it  is  either  the  case  that  three  are  positive  and  one  is  negative,  or  one  is  positive  and 
three  are  negative.  The  cases  are  generically  pictured  in  Figure  4.3. 

It  is  possible  to  determine  whether  pg{w )  has  four  real  roots,  or  two  real  and 
two  complex  roots,  by  hireling  the  sign  of  the  discriminant  of  pg{w )  (positive  in  the 
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(a)  Three  positive  (b)  Three  negative  (c)  One  positive,  one  negative, 

two  complex  conjugates 

Figure  4.3:  The  case  g  >  0.  Contours  of  constant  imaginary  part  equal  to  the 
imaginary  part  of  a  critical  point,  shown  for  x  E  ( — 7r,  7t]  . 

former  case,  negative  in  the  latter). 


4. 5. 3. 2.1  Constant  imaginary  part.  The  deformed  integration  path  should  be 
chosen  so  that  Irnag  (h(z))  is  constant.  We  have 

Irnag  (h(x  +  iy))  =  sin(x)  (2  cos(x)  [(<7/?2)e— 2y  —  (g(3-2)e2y] 

+  [(aq  +  gf3i)e~y  -  (am i  +  g(3-i)ey] )  (4.61) 

First,  consider  a  real  root  wq  of  pg(w).  Since  wq  is  real,  h(zo)  is  also  real,  so 

we  must  find  the  contours  where  Imag  (h(x  +  iy))  =  0.  This  equality  holds  when 

sin  (a;)  =  0,  or  when  cos(x)  =  —  ^ai^g,^)e-2t^gg+^2y^eV ,  which  approaches  0  as  y  — > 

±cx).  Therefore,  for  all  n  E  Z,  the  contours  x  =  mi  satisfy  Irnag  (h(z))  =  0,  and  the 

lines  x  =  |  +  nn  are  asymptotes  of  the  remaining  contours  with  Irnag  ( h(z ))  =  0.1 

1It  is  not  possible  to  choose,  e.g.,  x  =  0  as  the  deformed  contour,  because  the  deformation  would 
have  to  pass  through  regions  of  x- values  where  the  integral  does  not  converge. 
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Now,  consider  a  possible  complex  root  w0  of  pg(w).  We  must  find  the  contours 
where  Imag  (h(x  +  iy))  =  Imag  (h(zo)),  a  nonzero  constant.  Since  these  may  not 
cross  the  contours  where  Imag(h(z))  =  0,  they  are  restricted  to  regions  where  x  G 


(rm,  ( n  +  l)7r),  so  y  — »  ±cx).  Since  the  term  (g/3 2)e  2y  —  (gj3_2)e2y  dominates  the 


growth  in  y,  we  must  have  either  sin(x)  — >  0  or  cos(a;) 


(cvi+g/3i)e  y+(a-i+ff/3-i)ej/ 
2((gh2)e-2«-(g/3_2)e2i') 


fast  enough  to  keep  Imag  (h(z))  from  growing.  Once  again,  we  hnd  that  the  contours 


are  asymptotic  to  x  —  rm  and  x  =  |  +  nir  (n  G  Z). 


4. 5. 3. 2. 2  Steepest  descent  paths.  We  will  integrate  over  contours  where  the 
real  part  Real  ( h(z ))  reaches  a  maximum  at  one  of  the  critical  points.  We  have 

Real  (h(x  +  iy))  =  cos(2a;)  [(y/32)e_22/  +  (g/3-2)e2y] 

+  cos(x)  [(«!  +  gf3i)e~y  +  (a_ i  +  y/3_ i)ey]  +  (a0  +  g(30).  (4.62) 

Since  e±22/  dominate  the  growth  of  Real(h(z))  as  y  — >  ±oo,  and  /3±2  <  0,  we  have 
Real  (h(z))  — >  +oo  along  the  lines  x  —  |  +  rm  ( n  G  Z),  and  Real  (h(z))  — >  —  oo  along 
the  lines  x  =  rm  (n  G  Z).  We  now  have  three  cases  to  consider. 

In  the  first  case,  pg(w)  has  three  positive  roots  and  one  negative  root.  Let  w o  be 
the  middle  positive  root,  and  z0  =  —  iln(w0).  The  necessary  contour  in  Figure  4.3(a) 
joins  Zo  to  the  negative  root  and  its  reflection.  By  summing  equation  (4.54)  over  all 
adjacent  regions  of  width  27 r,  we  have 

T  1  1  ~1/2 

uj(t)  ~  U (zo  +  2nn )  '  ^  ~ 2/?//(^o  +  2nn )  '  wo°  ■  etHzo+2nn)t~1/2  (4.63a) 

nGZ 

_  r  i  i ~1/2 

=  y/n  --h"(z0)  ■  WQj  ■  e^^iT1/2  ■  '^'U(i  ln(w0)  +  2n7r),  (4.63b) 

^  nGZ 
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and,  because  U(—z )  —  U(z), 


U (i  In (w0)  +  2nir)  =  U(i  ln(w0))  +  y^  ( U (2m:  +  i  ln(w0))  +  U (—2m:  +  i  ln(w0))) 

nEZ  n= 1 

(4.64a) 

OO 

=  U(i  ln(w0))  +  2  Real  (7/ (2m:  +  i  ln(w0)))  G  M.  (4.64b) 

n=  1 

In  the  second  case,  pg(w )  has  one  positive  root  and  three  negative  roots.  Let 
Wq  be  the  middle  negative  root,  and  Zo  =  7T  —  In  | vjq \ .  The  necessary  contour  in 
Figure  4.3(b)  joins  z0  and  its  reflection  to  the  positive  root.  As  above,  we  have 
_  r  i  i  -1/2 

Uj(t)  ~  ■  WqJ  ■  eth^zoh~1/2  ■  ^  U(i  In  |iu0|  +  (2n  +  l)7r),  (4.65) 

^  nEZ 

and 

OO 

yy  U(i  In  |iy0|  +  (2 n  +  l)7r)  =  '^j(U((2n  +  1)7T  +  %  In  |iu0|) 

nSZ  ii=0 

+  U(—(2n  +  l)7r  +  i  In  |iu0|))  (4.66a) 

OO 

=  2  yy  Real  ( U((2n  +  l)7r  +  i  In  |w0|))  el.  (4.66b) 

n= 0 

In  the  third  case,  pg(w)  has  one  positive  root,  one  negative  root,  and  two  com¬ 
plex  conjugate  roots.  Let  w0  =  elZo  be  the  root  with  positive  imaginary  part,  and 
Wo  =  e the  root  with  negative  imaginary  part.  In  Figure  4.3(c),  the  necessary 
contour  for  zq  proceeds  from  the  asymptote  x  =  —  n  to  the  asympote  x  —  0,  while 
the  contour  for  —  zq  proceeds  from  the  asymptote  x  =  0  to  the  asympote  x  =  +7T. 
We  have 

T  1  1  ~1/2 

Uj(t)  ~  yy  U(zo  +  2m:)  ■  \[t:  —-h”(zo  +  2m:)  ■  Wq- 3  •  etft(2o+2^)^-i/2 

nEZ  L 
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r  i  i  ™1/2 

+  U(—zp  —  2to7t)  •  V7T  — -h"  (—zq  —  2mr) 

nEZ 

■  Wo™-'  •  eth(~*° -2n7r)^-1/2 

/  r  i  i™1/2 

=  v^t™1/2  ^  W(*o  +  2mr)  •  --/i"(2o)  V'e^2°) 

neZ  '  L 

_  f  1 _ .-1-1/2  _ \ 

+  U(z0  +  2nir)  •  ~-h"{z0 )  WJettW  j 

=  2v/vrt~1^2  ^  Real 

nEZ 

using  the  fact  that  h(—z)  =  h(z)  as  well. 


W(z0  +  2mr)-  -- h"(zQ ) 


eth(zo) 


(4.67a) 

(4.67b) 

(4.67c) 


4. 5. 3. 3  Relationship  between  the  solutions  (g  <C  1) 

Assume  that  a  root  w  of  pg(w)  can  be  expanded  as 

w  =  w0  +  gawi  +  g2aw2  H -  (4.68) 

for  some  a/0. 

First,  assume  that  Wq  =  0,  and  w  ~  gaW\.  Then 


g“‘+l{1hw‘1)+g3‘‘{a1w\)+  g°(-a-1w1)+  gH-ZjU)  ~  0-  (4.69) 

The  balance  between  (T)  and  (2)  yields  a  =  —  1,  and  thus  w\  =  —fh-  The  balance 
between  (3)  and  (4)  yields  a  —  1,  and  thus  W\  =  —  So  two  of  the  roots  are 


«i  _i  j  2(3-2 

w - —g  and  w - g. 

2  (32  o;_i 


(4.70) 


Now,  assume  that  w0  ^  0.  Then 


5f1(2/32Wg  -  2(3-2)  +  ^(opwo  -  a_iw0)  +  ga(3a1wlwl  -  a.- iiui)  ~  0.  (4.71) 
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The  O 


(1)  term  yields  Wq  =  The  balance  between  (l)  and  (5)  yields  a  —  1, 


and  thus  wi  =  —2  f2W°2  13  2  =  '3  2Ql,  32a  1 .  So  the  other  two  roots 

*  ociiuiq— a_i 


are 


w 


a_i 

CTi 


and 


w  ~ 


G!_l 

ai 


+  0Wi- 


(4.72) 


The  (9  (1)  roots  correspond  to  the  roots  in  the  deformed  contours  we  chose 
for  the  integration  —  they  are  either  the  two  complex  roots,  or  one  is  the  middle 
positive/negative  root,  since  the  remaining  roots  are  O  (g)  — >  0  and  O  (g^1)  —■ ►  ±oo. 


4. 5. 3. 3.1  Dependence  on  time  and  step  number.  The  step-step  interaction 
case  can  now  be  considered  as  a  regular  perturbation  of  the  no-interaction  case,  for 
steps  in  the  bulk  of  the  crystal. 

When  and  a_i  have  equal  signs,  which  is  possible  for  some  intermediate 
values  of  u0,  we  have  equation  (4.58c): 

Uj(t)  ~  v/^cqcmO-174  •  •  e^+v^fr1/2 

■  2  Real  (  U  (  (2 n  +  l)7r  —  i  In  ^  '^1 


n= 0 


The  time  dependence  is  proportional  to  e*(v/i°T+V/W-d)  which  is  exponential 

growth  with  a  polynomial  decay  term.  The  dependence  on  j  is  proportional  to 
(—  ,  which  is  growth  or  decay  depending  on  the  sign  of  j.  Note  that  the 

sign  of  Uj(t)  alternates  with  odd  and  even  values  of  j. 

When  q:i  and  a_i  have  opposite  signs,  which  holds  for  very  strong  and  very 
weak  electric  fields,  we  have  equation  (4.60c): 


u 


;(t)  ~  2v/vrt  ^  Real  j  U(zq  +  2nn)  ■ 


nEZ 


—2h'\z0) 


-1/2 


Jjh(z  o) 


w  be 
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First  consider  the  time  dependence.  Recall  from  equation  (4.59)  that  Real  (h(z0))  = 
Icq |  —  |a_i|.  We  found  in  Section  4.4.3  that  cq  decays  to  0  faster  than  a_i  when 


the  electric  held  is  very  strong,  implying  that  Real  (/r(20))  <  0.  Now  consider  the 


dependence  on  step  number.  Recall  that  w0  =  i 


so  again  we  have  growth  or 


decay  depending  on  the  sign  of  j. 
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4.6  Discussion 


We  have  analyzed  the  motion  of  crystal  steps  as  a  discrete  system  of  coupled  linear 
ODEs.  Much  of  the  literature  has  focused  on  continuum  approximations  of  surface 
height,  which  can  provide  reliable  indications  of  the  stability  of  uniform  step  trains; 
step  bunching  is  one  form  of  instability.  However,  these  continuum  approximations 
miss  a  time  dependence  term  that  appears  in  the  discrete  analysis,  and  they  do  not 
give  a  full  picture  of  the  dependence  on  step  number. 

Recall  our  measure  of  the  deviation  of  terrace  width  from  the  average  value, 
Uj(t)  =  xj{t)  _  i  We  found  that  the  time  dependence  of  Uj{t)  is  proportional 

to  eth^z°H~1/'2 .  The  f~1//2  factor  does  not  appear  in  continuum  approximations,  and  is 
the  fastest  decay  we  can  expect  from  a  first-order  asymptotic  analysis  of  the  discrete 
model.  The  h(zo)  exponent  depends  on  the  physical  parameters  of  the  problem. 
Stability  is  indicated  when  Uj(t)  decays  with  time,  while  instabilities  such  as  step 
bunching  are  possible  when  Uj(t)  grows  with  time. 

We  also  found  that  the  step  number  dependence  of  Uj(t)  is  proportional  to 
(e*2:o)"/.  Therefore,  Uj(t)  grows  or  decays  in  j  depending  on  the  sign  of  j.  We  have 
identified  parameter  regions  where  the  base  etz°  is  a  negative  real  number,  indicating 
alternating  signs  of  the  terrace  width  deviation.  Such  alternation  could  lead  to  the 
“double  steps”  that  have  been  observed  in  silicon  [60]  and  ruthenium  [64], 
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Appendix  A 


Addendum  to  Chapter  3 
A.l  Standard  ports 


In  2003,  many  peer-to-peer  applications  used  constant  port  numbers.  In  addition, 


many  client-server  applications  use  IANA-assigned  ports,  or  other  constant  ports. 


The  following  port  values  were  used  for  determining  the  “true”  classes  in  Section  3.3.1. 


Table  A.l:  Peer-to-peer  ports,  protocols,  and  their  usages.  Compiled  from 
various  online  sources. 


Port 

Protocol 

Application 

412 

TCP/UDP 

Direct  Connect 

1214 

TCP 

KaZaA 

1412 

TCP/UDP 

Direct  Connect 

2000 

TCP 

eDonkey  2000 

2004 

UDP 

eDonkey  2000 

2234 

TCP 

SoulSeek 

3415 

TCP 

KaZaA 

3531 

TCP/UDP 

KaZaA 

4329 

TCP 

iMesh 

4444 

TCP 

Napster 

4661 

TCP 

eDonkey  2000 

4662 

TCP 

eDonkey  2000 

4663 

TCP 

eDonkey  2000 

4665 

UDP 

eDonkey  2000 

4666 

UDP 

eDonkey  2000 

4670 

TCP 

eDonkey  2000 
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Table  A. 2:  Table  A.l  continued. 


Port 

Protocol 

Application 

5498 

TCP 

Hotline  Connect 

5499 

UDP 

Hotline  Connect 

5500 

TCP 

Hotline  Connect 

5501 

TCP 

Hotline  Connect 

5502 

TCP 

Hotline  Connect 

5503 

TCP 

Hotline  Connect 

5534 

TCP 

SoulSeek 

5555 

TCP 

Napster 

6257 

UDP 

WinMX 

6346 

TCP 

Gnutella 

6347 

UDP 

Gnutella 

6348 

TCP 

Gnutella 

6662 

TCP 

eDonkey  2000 

6666 

TCP 

Napster 

6699 

TCP 

WinMX 

6700 

TCP 

Napster 

6701 

TCP 

Napster 

6881 

TCP 

Bit  Torrent 

7674 

UDP 

Soribada 

7675 

TCP 

Soribada 

7676 

TCP 

Soribada 

7677 

TCP 

Soribada 

7788 

TCP/UDP 

BuddyShare 

8311 

TCP 

Scour 

8875 

TCP 

Napster 

8888 

TCP 

Napster 

8889 

TCP 

Napster 

22321 

TCP/UDP 

Soribada 

22322 

TCP 

Soribada 

41170 

TCP/UDP 

Blubster,  Piolet,  Manolito 
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Table  A. 3:  Client/server  ports,  protocols,  and  their  usages.  Compiled  from 
various  online  sources. 


Port 

Protocol 

Application 

20 

TCP 

FTP 

21 

TCP 

FTP 

22 

TCP 

SSH 

23 

TCP 

telnet 

25 

TCP 

SMTP 

53 

TCP/UDP 

DNS 

80 

TCP 

HTTP 

110 

TCP 

POP3 

123 

UDP 

network  time  protocol 

143 

TCP 

IMAP 

389 

TCP/UDP 

LDAP 

443 

TCP 

HTTPS 

993 

TCP 

IMAP  over  SSL 

1024 

TCP/UDP 

Windows  browsing 

1025 

TCP/UDP 

Windows  browsing 

1026 

TCP/UDP 

Windows  browsing 

1080 

TCP 

proxy  server 

1490 

TCP 

vocaltec  videoconferencing 

1755 

UDP 

windows  media  streaming 

1863 

TCP/UDP 

msn  messenger 

2703 

TCP/UDP 

SpamAssassin 

4099 

TCP 

AIM 

5050 

TCP 

Yahoo  IM 

5190 

TCP 

aim 

6277 

TCP/UDP 

DCC  anti-spam 

6665 

TCP 

IRC 

6667 

TCP 

IRC 

6670 

TCP 

vocaltec  videoconferencing 

7070 

TCP/UDP 

realaudio 

8000 

TCP 

Webserver 

8080 

TCP 

HTTP2 

22555 

UDP 

vocaltec  videoconferencing 

25793 

TCP 

vocaltec  videoconferencing 

32768 

TCP/UDP 

*nix  browsing 

32769 

TCP/UDP 

*nix  browsing 

32770 

TCP/UDP 

*nix  browsing 
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Appendix  B 


Addendum  to  Chapter  4 

B.l  Parameters  and  variables  of  the  crystal  problem 


Table  B.l:  Overview  of  parameters  and  variables. 


Symbol 

Dimensions 

Brief  description 

L 

length 

Typical  terrace  width 

a 

length 

Step  height 

0 

length2 

Atomic  volume 

Xj{t) 

length 

Position  of  jth  step  edge  at  time  t 

Uj(t ) 

0 

(xj+i(t)  —  Xj(t))/L ;  normalized  length  of  jth  terrace 

Uj(t) 

0 

Uj(t)  —  1;  deviation  from  uniform  terrace  width 

Co 

1/length 

Adatom  density  at  an  isolated  step 

Cj(x) 

1/length 

Adatom  density  at  position  x  on  the  jth  terrace 

Jj(x) 

1/time 

Adatom  current  (flux)  at  position  x  on  the  jth  terrace 

T 

time 

Desorption  time 

D 

length2 /time 

Terrace  diffusivity 

Vo 

length/time 

Down-step  velocity  clue  to  electric  held 

k+,  k- 

length/time 

Rate  coefficients  for  attachment-detachment 

9 

0 

Strength  of  step  interactions 

6 

0 

V 

0 

VDt/ (v0t  ) 
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